NASA TECHNICAL NOTE 




UJ 


JD 


LOAN COPY: RETUF = 
AFWL (DOGLJ 
KIRTLAND AFB, N. M. 


FITLOS: A FORTRAN PROGRAM FOR 
FITTING LOW-ORDER POLYNOMIAL SPLINES 
BY THE METHOD OF LEAST SQUARES 


by Patricia J. Smith 

Lewis Research Center 
Cleveland, Ohio 44135 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • AUGUST 1971 


TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 



1. Report No. 2. Government Accession No. 

NASA TN D-6401 

3. Recipient’s DIBSTIM 

4. Title and Subtitle FITLOS: A FORTRAN PROGRAM FOR FITTING 
LOW -ORDER POLYNOMIAL SPLINES BY THE METHOD OF 
LEAST SQUARES 

5. Report Date 

August 1971 

6. Performing Organization Code 

7. Author(s) 

Patricia J. Smith 

8. Performing Organization Report No. 

E -5292 

10. Work Unit No. 

129-04 

9. Performing Organization Name and Address 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio 44135 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Note 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D. C. 20546 

14. Sponsoring Agency Code 


15. Supplementary Notes 


16. Abstract 

FITLOS is a FORTRAN IV program to fit polynomial splines of degrees two and three. It 
combines some of the advantages of the method of least squares with the segmented curve 
of the theory of splines. FITLOS divides a set of data points into subsets and fits a poly- 
nomial of degree two or three on each subset by the method of least squares. The total curve 
is made smooth by making the polynomials on adjacent subsets and their first derivatives 
equal at the break point between the segments. For third-degree polynomials, the second 
derivatives are also made equal. These constraints are imposed by the method of Langrangian 
multipliers. This report describes the mathematical analysis of the least squares polynomial 
spline fit, gives complete documentation of program FITLOS, and serves as a user's guide 
for FITLOS. The program listing, flow charts, and example problems are included. 


17. Key Words {Suggested by Author(s) ) 


18. Distribution Statement 


Curve fitting 

Spline function 

Least squares polynomials 


Unclassified - unlimited 


19. Security Classif. (of this report} 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

93 

$3. 00 


For sale by the National Technical Information Service, Springfield, Virginia 22151 




CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION 1 

MATHEMATICAL DERIVATION 2 

Curve Fit 2 

Statistical Anaysis 8 

GENERAL DESCRIPTION OF PROGRAM 10 

HOW DATA ARE DIVIDED INTO SUBSETS 11 

INPUT DATA 13 

TYPICAL APPLICATIONS 16 

CONCLUDING REMARKS 18 

APPENDIXES 

A - PROOF THAT MATRICES X T WX AND C(X T WX) _1 C T 

HAVE INVERSES 19 

B - DETAILS OF SOLUTION OF EQUATION (5) 21 

C - PROGRAM LISTING AND FLOW CHART FOR FITLOS 28 

D - VARIABLES USED BY SEVERAL SUBROUTINES 35 

E - DESCRIPTION OF SUBROUTINES 36 

F - COMPUTER INPUT AND OUTPUT SHEETS FOR SAMPLE PROBLEM 1 ... 72 
G - COMPUTER INPUT AND OUTPUT SHEETS FOR SAMPLE PROBLEM 2 ... 81 

REFERENCES 91 


iii 


I 


FITLOS: A FORTRAN PROGRAM FOR FITTING LOW-ORDER POLYNOMIAL 
SPLINES BY THE METHOD OF LEAST SQUARES 
by Patricia J. Smith 
Lewis Research Center 

SUMMARY 

FITLOS is a FORTRAN IV program to fit polynomial splines of degrees two and 
three. It combines some of the advantages of the method of least squares with the seg- 
mented curve of the theory of splines. FITLOS divides a set of data points into subsets 
and fits a polynomial of degree two or three on each subset by the method of least 
squares. The total curve is made smooth by making the polynomials on adjacent sub- 
sets and their first derivatives equal at the break point between the segments of the 
curve. For third-degree polynomials, the second derivatives are also made equal. 
These constraints are imposed by the method of Lagrangian multipliers. 

FITLOS was written to complement other types of curve -fitting programs. This 
report describes the mathematical analysis of the least squares polynomial spline fit, 
gives complete documentation of the program FITLOS, and is intended to serve as a 
user’s guide for FITLOS. To augment this last purpose, the report includes examples 
of problems for which this type of curve -fit is useful. 


INTRODUCTION 

FITLOS was written to complement other curve -fitting programs. A new method 
of curve -fitting was needed that would combine some of the advantages of a least squares 
polynomial with the segmented curve of the theory of splines. Segmenting the curve 
gives it more freedom than a single polynomial over the entire range of the data, while 
fitting by the method of least squares smooths any small fluctuations in the data. 

The name ’’spline” is derived from the draftsman’s spline which is used to fair 
curves. Like the draftsman's spline, the spline function is smooth. DeBoor’s defini- 
tion of a spline function is used for this report (ref. 1). It is as follows: A function 
f(x) is a spline function of degree M with joints x^ < Xg < • • • <x fl if it has these two 
properties: 


L 


(1) In each of the intervals (-<*>, x^), J^x^, x 2 ), . . . , j^x n <x> ), f(x) is a polynomial 
of degree M. 

(2) The first M - 1 derivatives are continuous. 

In FITLOS, the continuity of the curve and its derivatives is imposed by the method of 
Lagrangian multipliers (ref. 2). 

The use of low -degree polynomials has two advantages. First, they have relatively 
few local maxima and minima. Second, they are easily differentiated and integrated. 
Second-degree polynomials have a third advantage; namely, their roots are easily 
found. Consequently, a FITLOS curve fit can be used readily for further applications. 

This report is intended to serve three purposes. First, it describes the details 
of the mathematical analysis of the least squares polynomial spline fit. Second, it pre- 
sents the program FITLOS, which makes this type of curve fit, and gives instructions 
for using the program. Third, it presents two problems for which the least squares 
polynomial spline fit is applicable and compares the results with fits made by other 
methods. 

Notation in the section MATHEMATICAL DERIVATION follows conventions in stand- 
ard mathematics textbooks. Involved proofs and mathematical details are given in the 
appendixes. 

To clarify the vocabulary, the word ’’order” refers to the sequence of points or 
numbers, while the word ’’degree” refers to the highest power of x in the polynomials. 
For example, the values x^ < x 2 < x^. are in order, while FITLOS fits polynomials 
of degree two or three. The difference between "subsets" and "segments” is a little 
more subtle. The set of data points is divided into subsets, while the fitted curve is 
divided into segments. However, the subsets of data correspond to the segments of 
the curve. 


MATHEMATICAL DERIVATION 
Curve Fit 

Consider a set of NX data points Z = |(x i; y.) | i = 1, 2, . . . , NX j where 
<x. For a weighted least squares polynomial fit of degree M, a matrix X can 
be defined: 
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Let W be the matrix of weights which has only diagonal elements, 


W = diag (w p w 2 , , w NX ) 


Let Y be the column vector 


Y - col (y p y 2 , , y^) 

where the y's have the same order as the x's in the matrix X. Let A be the column 

vector of undetermined coefficients. Then let Y be the column vector such that 
* 

Y = XA. For a weighted least squares fit, the scalar 

e - (Y - Y) 1 W(Y - Y) 

must be minimized with respect to each element of A. 

The weighted least squares polynomial spline fit can be described in a similar man- 
ner. First, however, a set of spline joints XM must be defined. Let XM be the 
set of x -values of the break points between the NS segments of the curve 
XM= |(xm) n ] n = 1, 2, . . . , NS - 1 j. Now set Z can be divided into NS subsets 
such that 


Z 1 = {( x i ! x i - x. < (xmjjJ 
Z 2 = {( x i> yp | ( xm )i 25 x j - ( xm ) 2 } 


Z NS " {/ X i’ y P | ^ xm ^NS-l - x i - x Nx} 
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For convenience, two sets of data point indices can be defined. Let F be the set of 
indices of the first data point in each subset F = {F„ }, where F is the smallest i 
such that (Xp y^) is an element of Z n . Similarly, let L be the set of indices of the last 
data point in each subset L = {L n }, where is the largest i such that (x^, y^) is 

an element of Z . From these definitions it can be seen that if any of the (xm) is an 

^ th ^ 

x-value of a data point, that data point is the last point in the n subset and the first 

th 

point in the (n + 1) L subset. However, the (xm) do not have to correspond to data points. 

When the data have been divided into subsets, a matrix X can be defined which is 
composed of submatrices such that 



for j = i, and X-^ is null for j ^ i. Matrix X has NS nonzero rectangular block sub - 
matrices on its diagonal and null submatrices elsewhere. The notation can be simplified 
a little at this point by dropping the second subscript on the submatrices of X since 
only diagonal elements are present. 

X = diag (X p X 2 , . . . X Ng ) 

Similarly, let Y be a column vector which is composed of NS subvectors Y. of the 
form 

Y i = 001 &F,' ^Fj+l’ ' ' • > 

Vector Y has the form 

Y = col (Yj, Y 2 , , Y NS ) 
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Let W be the matrix of weights which is composed of square submatrices of the 

form 


W ii = diag (w F w F v . . . , w L ) 
j 1 i 1 


for j = i, and W. . is null for i # i. 
J ’ i] J 

form 


Again dropping the second subscript, W has the 


W = diag (W 1 , W 2 , , W Ng ) 

Let A be a column vector of undetermined coefficients composed of NS subvectors of 
the form 


A j - co1 (a n , a i2> ■ • • > a i, M+l^ 

Vector A has the form 

A = col (A p A 2 , ... , A Ng ) 

Let Y be the column vector defined by the matrix product 

* 

Y = XA 


The scaler 


e = (Y* - Y) T W(Y* - Y) 


must be minimized with respect to each element of A, but subject to the constraints 
that the first M - 1 derivatives of Y* must be continuous at the break points between 
the segments of the curve. These constraints can be expressed in matrix form by de- 
fining the matrix C which is composed of submatrices 



H- 

2(xm). 
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for a quadratic fit, and 


c ir 


(xm). 

1 

0 


(xm)f 

2(xm) i 


3(xm)? 

6(xm). 


for a cubic fit for j = i. For j = i + 1, C, . = -C. . 1 . For other combinations of i 
and j, C^j is null. Again dropping the second subscript, C has the form 


C = 


C 1 -Cj 0 


0 C 2 -C 2 


C NS-1 C NS-1J 


The constraints take the form 


CA = 0 (1) 

The set of Langrangian multipliers can be introduced as a row vector A composed of 
NS - 1 subvectors of the form 

A i = ( A iP X i2’ ’ • ' ’ X iM^ 


Vector A has the form 


A = ( A p A 2 ; ' • • ’ A NS-1^ 

The scalar e then becomes 

e = (Y* - Y) T W(Y* - Y) + ACA 
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and must be minimized with respect to each element of A 
Substituting for Y* 

e = (XA - Y) T W(XA - Y) + ACA 

To minimize e with respect to a^, the derivative of e with respect to is set to 
zero; that is, 


0 = W(XA - Y) + (XA - Y) T W 9 i XA ~ Y ) + AC 


3a ii 3a ij 




1A 

8a i] 


Since 


S(XA - Y) 1 = 


8(XA - Y) 


9a iJ J 


W = W 1 , and a scalar is equal to its own transpose, we have 


S ( XA W(XA - Y) = (XA - Y) T W T d ^ X — — = (XA - Y) T WX 


9a ij 


Sa ii 


3a ii 


Therefore, 


2(XA - Y) WX + AC 


5A 

3a ij 


= 0 


(2) 


Equations (1) and (2) must be solved for A and A. Since 3A/9a^. # 0 for any a^, it 
must be true that 


2(XA - Y) T WX + AC = 2A T (X T WX) - 2(Y T WX) + AC = 0 


(3) 


rp T - 1 

Since the matrix X WX has an inverse, right multiplying by (X WX) , dividing by the 

T 

scalar 2 , and separating the unknowns A and A gives 


A T + — AC (X T WX)" 1 = ( Y T WX) (X T WX) " 1 


(4) 


L 



Ill 


iiiiaiiHw 


i 


in* mi ■ i i ■■ hi i 


I JIMI 


II I I I II I 


T 

The proof that X WX has an inverse is given in appendix A. 

Since CA = 0, (CA)^ = A T C T = 0. Right multiplying equation (4) by C T gives 


- AC(X T WX)" 1 C T = (Y T WX)(X T WX) 1 C T 
2 


Since the matrix C(X r ^WX)”^C^' 
gives 


has an inverse, right multiplying by 


C(X T WX) X C T 


1 A = (Y T WX)(X T WX)‘ 1 C T 

2 


C(X T WX) 1 C T 


The proof that C(X T WX) _1 
Substituting for 1/2 A 


T 

C has an inverse is also given in appendix A. 

T 

in equation (4) and solving for A gives 


A T = (Y T WX) l 1 - (X T WX)“ 1 C T 


C(X T WX) 1 c T 


-1 


C J (XWX) 


-1 


(5) 


The details of this matrix manipulation and a method of finding 
given in appendix B. 


C(X T WX) 1 C T 


is 


Statistical Analysis 


FITLOS makes a rudimentary statistical analysis of the curve -fit. It calculates the 
deviation and relative error between the given data and the fitted curve, the variance 
and the standard deviation, and Pearson's correlation coefficient. The formulas were 
taken from reference 3, but they are standard in any statistics textbook. The formulas 
are as follows; 

Deviation: 


d. 

i 



- Vi 


Relative error: 


e i 


d. 


l 
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Variance: 


where 


and 


°2 ~ 


NX 

Z < d i - 3 > 

i=l 


3\2 


NX 

E ^ 


d = 


NX 


F = Number of degrees of freedom 

= Number of points - Number of constraints 
- NX - M(NS - 1) 

Standard deviation: 

Correlation coefficient: 
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Since the number of constraints is large, and hence, the number of degrees of free- 
dom is small, the correlation coefficient can be deceptively small. For this reason, 
FITLOS also calculates the maximum possible correlation coefficient, which is r for 
y* = y^ for all i. The maximum r is equal to F/NX. 


GENERAL DESCRIPTION OF PROGRAM 

FITLOS was written in FORTRAN IV for the computer at the Lewis Research Center, 
which is an IBM 7094 11/7044 or 7040 Direct Couple computer under IBSYS version 13 
using ALTIO. 

Computer storage for the program as it is presented here, with 350 data points and 
10 spline joints, is around 20 000 locations. Since the Lewis computer has 32 000 loca- 
tions, the program could be expanded to fit more data points or to fit the curve in more 
segments. 

The program is written as a series of subroutines so the actual curve -fitting routine 
could be used as part of another program. The actual fit requires only three subroutines: 
one to divide the data into subsets, one to define the matrices in equation (5), and one to 
solve equation (5). 

In order to make the subroutines as flexible as possible, their arrays have variable 
dimensions. To conserve execution time, every subprogram with variably dimensioned 
arrays is called only once by its subroutine name. These calling vectors contain only 
the array names and the dimensions of the corresponding arrays in the main program 
FITLOS. Afterwards, the subroutines are called by entry names which do not disturb the 
size of the variably dimensioned arrays set up by the first call by the subroutine names. 

The main program FITLOS reads input data, calls the subroutines to make the fit 
(see tree diagram for hierarchy of subroutines, fig. 1), makes a statistical analysis of 
the fitted curve, and writes the output data. 


FITLOS 



— 

J 


Calls to set dimensions in subroutines 



j 

0RD, SEG, DEF, SLV, RFT, TRN, DXM, 


TRANSF j 

DSPL, DH, DL, and DNS 

1 



r~ 

|divxm 



Figure 1. - Tree diagram of subroutine calls. 
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FITLOS uses the following procedure to fit a curve: It reads input data, which are 
described in the section INPUT DATA. After the data are read, FITLOS checks logical 
variables TRANX and TRANY. If either is . TRUE. , subroutine TRANSF is called to 
make a log transformation on x or y. 

The next subroutine called is ORDER, which arranges the data points in order of 
ascending x. The next subroutine called is SEGMNT, which divides the data into sub- 
sets. SEGMNT is a monitoring routine which controls calls to small subroutines 
(DIVXM, DIVNS, FFLOW, FFHIGH, and SPESHL) which actually allot the data to subsets. 
There are four methods of dividing the data into subsets. These are described in the 
section HOW DATA ARE DIVIDED INTO SUBSETS. 

When the data have been divided into subsets, subroutine DEFMAT is called to de- 
fine the matrices (X^WX) - ^, Y^WX, and C. Then subroutine ASOLVE is called to 
perform the matrix manipulation involved in solving equation (5). 

FITLOS can check whether the curve was fit in more segments than were necessary. 

If the logical variable LREFIT is . TRUE. , subroutine REFIT is called to do this check- 
ing. When the refit checking is complete, FITLOS again interrogates logical variables 
TRANX and TRANY. If either is . TRUE. , subroutine BTRANS is called to transform 
the data back to its original form. Then the statistical analysis is made and the output 
data are written with descriptive labels and headings. 

If REFIT has indicated there were too many segments in the first fit, new spline 
joints are determined, subroutines DEFMAT and ASOLVE are called again, a new statis- 
tical analysis is made, and the new output data are printed. 

A listing of each subroutine, along with a flow chart and a description of its opera- 
tion, is provided in appendix E. Variable names and their limitations or special fea- 
tures are found in the, program listings. More details of how FITLOS works can be 
found in the section INPUT DATA. 

HOW DATA ARE DIVIDED INTO SUBSETS 

FITLOS provides four methods of dividing the data into subsets. The user deter- 
mines which method is used by proper setting of the input variables. 

The user has the option of selecting the spline joints, of selecting the number of 
segments, or of choosing one of the methods the program does automatically. 

If the user selects the spline joints, he must supply these data as part of the input. 
Then subroutine DIVXM searches the x -array to determine the index of the first and 
last point in each subset. 

If the user chooses the number of segments, subroutine DIVNS is called to divide the 
data as evenly as possible among the subsets. If the user does not specify the number 
of segments, the program will automatically choose the largest possible number of seg- 


11 


merits based on the number of data points and the degree of the polynomials. Again, sub- 
routine DIVNS is used to divide the data into subsets. 

If the user does not specify either the spline joints or the number of segments, 
subroutine SEGMNT checks the number of data points. If there are less than 3M + 1 
points, subroutine SPESHL is called to make a special division of the data into subsets. 
For M = 2, the division is as follows: 


Index of 

Index of 

Number of 

Spline joints 

first points 

last points 

subsets 


1 

3 or 4 

1 

x 3 or x 4 

1,3 

3, 5 

2 

X 3’ X 5 

1,3 

4,6 

2 

x 3 + x 4 ^ 

2 ’ 6 

For M = 3, SPESHL divides the data as follows: 


Index of 

Index of 

Number of 

Spline joints 

first points 

last points 

subsets 


1 

4, 5, or 6 

1 

x 4 ,x 5 , or x 6 

1,4 

4,7 

2 

X 4’ X 7 

1,4 

5,8 

2 

X 4 + X 5 x 
2 ’ 8 

1,5 

5,9 

2 

X 5’ X 9 


The final method of dividing the data into subsets is by force-fitting. The first 
M + 1 data points are used to determine a Lagrange interpolation polynomial. If the 
next point, the (M + 2)" 1 point, falls on the polynomial, it is accepted in the first sub- 
set. Then the next point is examined, and so on to the end of the set of data points. If 
a point does not fall on the polynomial, a new subset is started with the next M + 1 
points. There are two subroutines to do a force-fit division of the data. FFLOW starts 
at the low end of the data set and FFHIGH starts at the high end. 
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INPUT DATA 


Input to FITLOS is by punched cards. The order of these cards, their formats, the 
variables they contain, and the use of these variables in controlling how the curve is fit 
are as follows: 


Card 

Format of card 

Variable 

Description 

1 

(12A6) 

TITLE 

Alphanumeric identification of the 
data. The title must be confined to 
columns 1 to 72 of one card. 

2 

(5I3,4L3,F12. 6) 

M 

Degree of the polynomial. M must 
be 2 or 3. 



NX 

Number of (x, y, w) data points. 



NS 

Number of segments if the user 
selects the number of segments. 

NS * 0 means the curve will be fit 
in NS segments. If NS = 0, the pro- 
gram will select the largest possible 
number of segments. 



NB 

Number of spline joints. NB + 0 
indicates the user has selected the 
spline joints and these data will be 
read as part of the input data. 

NB = 0 means the program will set 
NB = NS - 1. 



NF 

Numerical variable which indicates 


whether force -fitting is used to di- 
vide the data into subsets. If N'F <0, 
the data are divided into subsets by 
force-fitting starting at the low end 
of the data. If NF > 0, the data are 
divided by force -fitting starting at 
the high end of the data. If NF = 0, 
the program divides the data as 
evenly as possible among the maxi- 
mum possible number of subsets. 
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Card 


Format of card 


Variable 


Description 


2 (5I3,4L3,F12. 6) LREFIT Logical variable which indicates if 

FITLOS should check whether the 
curve was fit in more segments than 
were necessary. LREFIT = .TRUE, 
means a check should be made. 
LREFIT = . FALSE, means no check 
should be made. The write-up of 
subroutine REFIT gives details of 
how the check is made. 

TRANX Logical variable which indicates if a 
log transformation should be made 
on x and (xm). TRANX = .TRUE, 
means the transformation should be 
made. TRANX = . FALSE, means 
the transformation should not be 
made. 

TRANY Logical variable which indicates if a 
log transformation should be made on 
y and y*. TRANY = .TRUE, means 
the transformation should be made. 
TRANY = .FALSE, means the trans- 
formation should not be made. 

NPUNCH Logical variable which indicates if 

cards containing the coefficients 
should be punched. NPUNCH = 

. TRUE, means no cards should be 
punched. NPUNCH = . FALSE, 
means cards should be punched with 
all the coefficients for one segment 
on a card. 

TOL Tolerance acceptable for refit 

checking or for force -fitting. De- 
tails of how TOL is used are found 
in the descriptions of subroutine 
FFLOW. 
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Card 

Format of card 

Variable 

Description 

3 

(12A6) 

FMT 

Variable format for reading (x,y,w) 
data points. 

4 to 3 + n (n = 

(FMT) 

X 

Independent variable array. 

number of data 


Y 

Dependent variable array. 

cards) 


W 

Array of weights. Since FITLOS 
makes a weighted least squares fit, 
each point must have a weight. How- 
ever, if all the weights are zero, 
FITLOS will make all the weights 1. 

The (x, y,w) data are read in the order (x^, y^, w 

l), (x 2 , y 2 > w 2^’ • • • » ( x NX> y NX’ 

Wnx ). If NB y 0, 

the following data are needed: 


4 + n 

(12A6) 

FMTM 

Variable format for reading the spline 
joints selected by the user. 

5+n to 3+n+m 
(m = number of 
spline joint 
cards) 

(FMTM) 

XM 

Array of spline joints. If NB = 0, 
FMTM and XM are not needed. 

4 + n + m 

(13) 

KASES 

The number of additional fits to be 
made with the current (x, y,w) data. 


KASES is originally set to zero by 
FITLOS so a title card and (x, y, w) 
data are read in. If the KASES card 
contains zero or is blank, FITLOS 
will transfer to read a new title card 
and new (x, y, w) data. If KASES y 0, 
FITLOS will transfer to read a new 
card 2. KASES is reduced by 1 each 
time a new card 2 is read in until 
KASES finally becomes 0. 

Variables NS and NB are not independent. FITLOS interrogates NB to determine if 
more input cards should be read. Subroutine SEGMNT interrogates NB first. If 
NB y 0, NS is set equal to NB + 1, and the division into subsets is based on NB and the 
chosen spline joints. If NB = 0, SEGMNT interrogates NS. If NS y 0, NB is set equal 
to NS - 1, and the division into subsets is based on NS. If both NB and NS are zero, 
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SEGMNT sets NS equal to the maximum possible number of subsets and then interrogates 
NF. More specific details of how these input variables are used can be found in the de- 
scriptions of the individual subroutines. 


TYPICAL APPLICATIONS 

One typical application for FITLOS is fitting- experimental data. An example of 
this is the calibration of a multiplier phototube -capacitor, where the independent varia- 
ble is time and the dependent variable is digitizer counts. Data from several different 
light sources are translated until the curves coincide as nearly as possible. Since the 
curves do not coincide exactly, there are small fluctuations in the data. For such a 
calibration to be useful, these fluctuations must be eliminated. 

Obviously, any least squares fit would do that. However, fitting this curve with a 
single least squares polynomial of degree one, two, or three did not give satisfactory 
results. Figures 2 to 4 (appendix F) show the relatively large deviations between the data 
points and the fitted curve. The curve was then fit using FITLOS with three polynomials 
of degree two. The deviations between the data points and the fitted curve are sufficiently 
small, as figure 5 shows. 

The curves in figures 2 to 5 (appendix F) are plotted on a log -log scale to emphasize 
these deviations. The plots of the deviations are made on a semi -log scale because they 
are both positive and negative. The computer listings from which these plots were made 
and the computer input sheet for the FITLOS fit can be found in appendix F. 

Another application for FITLOS is approximating a curve to obtain further informa- 
tion about it, such as the derivative and the definite integral. The source of the data 
points is immaterial. They could be experimental data points or they could be generated 
from some complicated function. The points for this example were generated from the 
function 


f(x) = x sin x - 1 

This function was chosen for the example because it is not a polynomial and yet it is 
simple enough to be differentiated and integrated analytically for comparison with the 
results from FITLOS. The derivative of f(x) is 

f'(x) = x cos x + sin x 


and the definite integral is 


/ 


X 


f 


f(x)dx = (sin x - x cos x - x) 


For finding the derivative and the definite integral using a FITLOS curve, the third - 
degree polynomials yield smoother curves. The derivative is 

2 

y*’ = a 2 + 2a g x + 3a 4 x 

The integral is a little more complicated because each segment of the curve must be in- 
tegrated separately. Consequently, the definite integral takes the form 


x f Ax m). UN ,(xm) n 

J y* dx = J y* dx+ J y n ^ + J 

x x„ n=i+l (xm)„ 1 (xrr 


•( xm ) 


r * dx 


y i+N 


n=i+l (xm) n _ 1 


(xm). 


i+N 


where i is the number of the first spline joint such that (xm). -l <x o< (xm)p and N 
is the number of additional segments such that i + N s NS and (xm)^ < x^. Tables I 


to III (appendix G) compare the FITLOS curve y* with f(x), y*» with f'(x), and / y* dx 

/* x f X o 

/ f(x) dx. Figures 6 to 8 (appendix G) are plots of the data in these three tables. 


with 


To determine the roots of this curve, the curve should be fitted with second-degree 
polynomials. The roots of y* can be found by the quadratic formula. The roots of 
f(x) can be found numerically (by the Newton -Raphson method) for comparison. The 
following table compares the Newton -Raphson roots with the FITLOS roots: 


Newton-Raphson root 

FITLOS root 

Deviation 

1. 1141571 

1. 1143261 

-0. 0001690 

2. 7726047 

2. 7714741 

0. 0011306 


The computer listings and the input sheet for FITLOS for this example are pre- 
sented in appendix G. 
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Another application for FITLOS, one that is shared by all curve -fitting schemes, is 
generating points for mechanical plotting. Automatic plotting devices such as the 
Calcomp Plotter or the DD80 Microfilm Plotter require a method of generating points 
close together. Figures 2 to 5 illustrate this application, since these plots were done 
on the Calcomp Plotter at the Lewis Research Center. Figures 6 to 8 were done on the 
DD80 Microfilm Plotter. 


CONCLUDING REMARKS 

This report has described the mathematical analysis of the least squares polynomial 
spline method of curve fitting; has presented the FORTRAN program FITLOS, which 
makes this type of curve fit; and is intended to serve as a user’s guide for FITLOS. The 
sample problems included show problems for which this type of curve fit is useful. 

They also show how the curve fit may be used for further applications such as integration, 
differentiation, root finding, and plotting. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 23, 1971, 

129-04. 
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APPENDIX A 


PROOF THAT MATRICES X T WX AND C(X T WX) _ 1 C T HAVE INVERSES 

Matrix X^WX 

Let X and W be the matrices in the main-text section Curve Fit. Since X is 

T 

block diagonal and W is diagonal, the product matrix X WX is block diagonal with 
diagonal blocks of the form 



Then (X T WX)i becomes (X T WX). - X: F U T UX.. If P is defined as P = UX., then 

rp _T 111 1 

(X A WX) i = P T P. 

Since the leading principal minor of X i is Vandermonde of order M + 1 and since 

none of the Xj are equal, by the definition of the spline function, Xj has rank M + 1. 

Since premultiplying X^ by the nonsingular matrix U does not change the rank, the 

product P has rank M + 1, by theorem 5. 6. 3 of reference 4. The matrix 
T T 

(X WX)i = P P then has rank M + 1, by theorem 5. 5. 4 of reference 4. Therefore, 
since (X T WX). has dimension (M + 1) x (M + 1) and has rank M + 1, it is nonsingular. 
Consequently, the inverse (XWX) i exists. 
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m T — 1 

Since (X WX). is defined for all i, all the submatrices (X WX). exist and, hence, 
i T -1 1 

the entire inverse (X WX) exists. 

Matrix C(X T WX) _1 C T 

T -1 

Let (X WX) be the inverse matrix found in the preceding section. Let C be 
the matrix of constraints defined in the main-text section Curve Fit. Since has 
rank M and since there are NS - 1 rows of blocks in C, the rank of C is (M)(NS - 1). 

«T 

It was shown that X WX is positive definite since it can be decomposed into the 
form 


(X T WX) = P T P 

T 1 T -1 

Consequently, its inverse (X WX)“ is also positive definite. Therefore, (X WX) 

_ p 

possesses a positive definite square root Q (see pp. 92 to 93 of ref. 5), and (X WX) 
can be written as 


(X T WX) _1 = QQ T 
T -1 

where Q has the same rank as (X WX) , which is (M + 1)(NS). Therefore, the matrix 
C(X T WX) _1 C T can be written as 

C(X T WX) _1 C T = CQQ T C T 

or as 

C(X T WX) _1 C T = PP T 

where P = CQ. Since postmultiplying C by the nonsingular matrix Q leaves the rank 
of the product unchanged, P has the same rank as C, which is (M)(NS - 1), by theorem 
5. 6. 3 of reference 4. The matrix C(X T WX) _1 C T = PP T then has the rank (M)(NS - 1) 
by theorem 5. 5. 4 of reference 4. 

Since CCX^WX)”^^ also has dimension (M)(NS - 1) x (M)(NS - 1), it is nonsingular. 
Therefore, its inverse [c(x T wx ) _1 C T ] exists. 
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APPENDIX B 


DETAILS OF SOLUTION OF EQUATION (5) 

The solution of equation (5) requires some rather involved matrix manipulation. 

The calculation of the matrix C(X^WX) and its inverse is particularly complicated. 

Let us define the matrix B to be 

B = C(X T WX) _1 C T (6) 

and its inverse to be B = D. From the definition of the partitioned matrices C and 
(X WX) , it can be seen that B is composed of submatrices of the form 

B ij = c i [(x T wx): 1 + (x T wx)- + 1 1 ]c: r 

for j = i and i = 1, . . . , NS - 1; 

B.. = -C i (X T WX)“ 1 cJ r 

for j = i - 1 with i = 2, . . . , NS - 1 and for j = i + 1 with i = 1, . . . , NS - 2; 


for other combinations of i and j. Since there are at most only three nonzero subma- 
trices in each row of B, these can be redefined as follows: 


B il = -C i (X T WX): 1 cT_ 1 for 

i= 2, . . . 

, NS - 1 

B i2 = C.f^wx): 1 + (xTwxj^Jc^ 

for i = 1, 

. . . , NS - 1 

B i3 = -Cj(X T WX)r^ lC T j for 

i= 1, • • • 

, NS - 2 


Consequently, matrix B takes the form 
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B 21 B 22 B 23 0 
0 B 31 B 32 B 33 0 


B = 


0 ► 

B NS-3, 1 B NS-3, 1 

B NS-3, 3 

0 

0 — -► 

0 B NS-2, 1 

B NS-2, 2 

B NS-2, 3 

0 

0 

B NS-1,1 

B NS-1, 2 


Matrix D has the form 


D = 


D 


11 


D 


NS-1, 1 


D 1,NS-1 


d ns-i,ns-i 


From the definition of the inverse of a matrix, BD = I, it can be seen that the product 
matrix BD has the form 


( B 12 D 11 + B 13 D 2l) 

( B 21 D 11 + B 22 D 21 + B 23 D 31^ 
( B 31 D 21 + B 32 D 31 + B 33 D 41^ 


BD = 


( B 12 D 12 + B 13 D 22^ 


(B NS-2, l D NS-3, l + B NS-2, 2 D NS-2, 1 +B NS-2, 3 D NS-1, 1^ 
^ B NS-1, l D NS-2, 1 + B NS-1, 2°NS-1, 1^ 
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If I is partitioned into submatrices of the form 6^, where 6^ is the identity 
matrix for j = i and 6^ is null for j ± i, then each column k of the product matrix 
BD becomes a series of simultaneous linear matrix equations of the form 


B 12 D lk + B 13 D 2k = 6 lk 

(7) 

+ B 22 D 2k + B 23 D 3k = 6 2k 

(8) 

+ B 32 D 3k + B 33 D 4k = 6 3k 

(9) 


B NS-2, l D NS-3, k + B NS-2, 2 D NS-2, k + B NS-2, 3 D NS-1, k ~ 6 NS-2, k 

B NS-1, l D NS-2, k + B NS-1, 2 D NS-1, k = 6 NS-1, k ( 10 ) 


,-l 


The solution of this system is begun by left multiplying equation (7) by (The proof 

that B^ has an inverse is the same as the proof for the existence of (jCCX^WX)’"'^^! 
found in appendix A. ) Solving equation (7) for gives 


D lk “ B 12 6 lk " B 12 B 13 D 2k “ B 12 ( 6 lk B 13 D 2k 


Substituting for in equation (8) and solving for gives 


°2k “ ( B 22 " B 21 B 12 


2 B l 3 ) _1 ( 


6 2k " B 21 B 12 6 lk " B 23 D 3k) 


Similarly, substitution of into equation (9) will give a similar solution for D^. 

This process can be repeated for the entire set of equations. 

However, the matrix algebra can be simplified by defining two auxiliary matrices 
E and A. Let E^ = B^ and let = B^g _ g for 1 + 1. Let 

and let A^ = 6^ - B^jE”^ f° r 1*1- Then the solution of the first NS - 2 

equations can be written as 
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( 11 ) 


D £k -E £ ( A l ~ B l,3 D l+l,k ) 

Substituting for j ^ into equation (10) and solving for D^ s ^ k gives 


D NS-l,k~ E NS-1 A NS-1 


Since all the E^ and can be found in terms of known quantities, D^g ^ k can be 
determined uniquely. Similarly, all the D^ k can be found by substitution into equa- 
tion (11) for l = NS - 2, . . . , 1. 

This scheme is easily programmed. Since the B.j have dimensions 3X3 for a 
quadratic and 4X4 for a cubic, the matrix E is either a 3 x 3 or a 4 x 4. Consequently, 
the largest matrix that must be inverted by numerical methods is a 4 x 4. This inver- 
sion can be done with good accuracy by any standard numerical matrix inversion tech- 
nique. Even though this scheme involves many arithmetic operations, the round-off 
error in the final answers becomes apparent only in the 12th or 13th significant figure. 

Once the matrix D is determined, equation (5) becomes 


a t = (yWi - (X T WX)" 1 C T DC 


(X T WX)" 1 


T 

The next problem is to form the matrix product T = C DC. This is somewhat compli- 
cated because the submatrices of T take special forms depending on their location in 
the matrix. These forms are readily seen from the definition of the partitioned matrices 
C and D. The forms are summarized as follows: 


T. = s(cTd. C ^ 

]« V J jn n) 

for" the corner elements. The scalar s equals +1 if n=j and -1 if n * j. 

T. = sC T (D. C - D. ,C ,) 

]n ]' jn n jn-1 n-V 

for the noncorner elements of the top and bottom rows. The scalar s equals +1 for 
j = 1 and -1 for j = NS. 


T. = s(c T D. - C T D. 1 )c 
jn \ ] jn jn j-l,ny n 
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for the noncorner elements of the first and last columns. The scalar s equals +1 for 
n = 1 and -1 for n = NS. 


T i,n= ( c i'-l D j-l,n-X - C f D J-l,n-lK-l + ( C ?' D j. 



j,n-l 



for the elements in the "middle” of the matrix. 

Once T is defined, the solution of equation (5) is quite straightforward. Since 
the matrices are partitioned, the matrix multiplication can be done in several steps 
which can be programmed easily. These steps are as follows: Since equation (5) 
becomes 

a t = (Y T WX)[i - (x t wx) _1 t (x t wx) _1 

the first step is to carry out the multiplication by Y WX. Equation (5) becomes 

A T = (Y T WX) - (Y T WX)(X T WX)" 1 T (X T WX) _1 
Letting the product (Y^ WX)(X T WX) _1 define the vector VV, equation (5) becomes 

A T = [(Y T WX) - (VV)t] (X T WX)~ 1 

Letting the vector V be defined as (Y T WX) - (VV)T, equation (5) becomes 

a t = V(X T WX) _1 

Writing out each of these steps in terms of the partitioned matrices and vectors gives 

(X T WX)^ 1 

A T = (A? A 5 S ) = (v v V NS ) 

0 

_ 

= [v^xTwX)- 1 , . . . , V NS (X T WX)- 1 g 

T 

Consequently each subvector of A can be determined independently of the others. 

ri rp 

The n subvector of A can be written as 
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A n = V^WX)' 1 


Writing out V in terms of its definition gives 


V = (Y i WX) - (VV)T 


( V 1 V Ns)= [( yTwx > 1 - < yTwX )n S ] - ( VV 1> • • ' > ^Ns) 


T n • • • T 1,NS 


l_ T NS,l • • • T NS,NSj 


NS NS 

(Y TwX )l - Z VV j T il- ' ■ ■ > < yTwX >NS - Z VV j T iNS 
j=l j=l 


Therefore, V becomes 
’ n 


... NS 

V„= < Yl WX) n - £ V Vj T jn 

J=1 


Writing out VV in terms of its definition gives 


VV - (Y T WX)(X T WX) 1 


(Wj, . . . , vv NS ) = [(y t wx) 1? . . . , (Y T WX) 


NS 


(x T wx): 1 


(X T WX)^ S 


- [(Y T WX) 1 (X T WX ) ; 1 , .... (Y T WX) NS (X T WX)- 1 g ] 


Consequently, 


VV. = (Y T WX).(X T WX). 1 


T t -1 

Writing (Y WX)j, (X WX)^ , and T^ n in terms of their row and column elements, VW 
becomes 
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V V [(Y T WX) m , . . . , (Y T WX) j)Mtl ] 


( VV J,1- • • • - VV J,M + l) 


< xTwx >U, 1 • • • (x T wx)r 1 i^ M+ i 


L^w*))'; Wi, i • ■ • ( xTwx )[,m + i,m + iJ 


M+l 

E 

k-1 


(Y T WX) jik (X T WX)- 1 k l> . 


M+l 


■ . E ( YTwx )),k( xTwx >[,k, 


M+l 


te=l 


Similarly, V n becomes 


NS 


V n =(Y T WX)„- J2 ( Y \l vv j;M+1 ) 

j=l 


T- 


j, n, 1, 1 • 


T j,n,l,M+l 


T. 


- T, 


j,n, M+l, 1' ' ‘ j,n, M+l, M+l 


( V n,l > V n,M + l) - [< yTwX > M > ' ' ' ' (YTwX) n, 


M+l 


/ M+l NS 

~ / j ( 2 VV j,k T j, n,k, V ■ ’ X/ V j,k T j,n,k, 
\ k=l k=l 


M+l 


NS M+l 

(Y T WX) -E E VV i,k T j,n,k,l>- >< YTwX )n,M + l 

3=1 k=l 

NS M+l 

" 2 Z VV j,k T j,n,k,M+l 
j=l k=l 


Once the vector V has been formed, it is a simple matter to combine it with the sub- 
n ’ 

matrix (X^WX)”'*' to get A^. 
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APPENDIX C 


PROGRAM LISTING AND FLOW CHART FOR FITLOS 


IBFTC FITLOS 

INPUT VARIABLES 


TITLE - HOLLERITH IDENTIFICATION OF PROBLEM 

M - DEGREE OF POLYNOMIAL 

NX - NUMBER OF DATA POINTS 

NS - NUM8ER OF SEGMENTS CHOSEN BY USER 

NB - NUMBER OF SPLINE JOINTS (INTERIOR) IF USER CHOOSES THEM 
NF - NUMERICAL SIGNAL TO DETERMINE WHICH AUTOMATIC METHOD OF 

DIVIDING THE DATA INTO SUBSETS - USED ONLY IF BOTH NS 10 

AND NB ARE ZERO 11 

L REF I T - LOGICAL VARIABLE 12 

IF LREFIT IS TRUE, PROGRAM WILL CHECK FOR DUPLICATION 13 

OF COEFFICIENTS 16 

IF LREFIT IS FALSE, NO CHECK WILL BE MADE 15 

TRANX - LOGICAL VARIABLE 16 

IF TRANX IS TRUE, A LOG(IO) TRANSFORMATION WILL BE 17 

MADE ON X AND XM 18 

IF TRANX IS FALSE, NO TRANSFORMATION WILL BE MADE 19 

TRANY - LOGICAL VARIABLE 20 

IF TRANY IS TRUE, A LOG(IO) TRANSFORMATION WILL BE 21 

MADE ON Y AND YC 22 

IF TRANY IS FALSE, NO TRANSFORMATION WILL BE MADE 23 

N PUNCH - LOGICAL VARIABLE 26 

IF NPUNCH IS TRUE, NO COEFFICIENT CARDS WILL BE PUNCHED 25 

IF NPUNCH IS FALSE, SEGMENT COEFFICIENTS WILL BE 26 

PUNCHED ON CAROS 27 

X - ARRAY OF INDEPENDENT VARIABLES 28 

Y - ARRAY OF DEPENDENT VARIABLES 29 

W - ARRAY OF WEIGHTS - MAY BE READ AS ALL ZEROS 30 

XM - ARRAY OF SPLINE JOINTS 31 

TOL - TOLERANCE FOR FORCE FITTING AND REFIT CHECKING 32 

KASES - NUMBER OF ADDITIONAL CASES USING SAME DATA 33 

36 

VARIABLES USED IN SUBROUTINE CALLS 35 

36 

37 

XX - ARRAY OF ORDERED INDEPENDENT VARIABLES 38 

YY - CORRESPONDING ARRAY OF DEPENDENT VARIABLES 39 

WW - CORRESPONDING ARRAY OF WEIGHTS 60 

NXX - NUMBER OF POINTS IN XX, YY, AND WW ARRAYS 61 

XM - ARRAY OF SPLINE JOINTS 62 

LLOW - ARRAY OF INDICES OF FIRST POINTS IN EACH SUBSET 63 

LHIGH - ARRAY OF INDICES OF LAST POINT IN EACH SUBSET 66 

XWX - MATRIX (X-TRANSPOSE«W»X!-INVERSE 65 

YWX - VECTOR ( Y-TRANSPOSE*W*X) 66 

C - MATRIX OF CONSTRAINTS 67 

A - VECTOR OF UNDETERMINED COEFFICIENTS 68 

YC - ARRAY OF DEPENDENT VARIABLES CALCULATED FROM EOUATION 69 

YC = XA 50 

51 

PROGRAM VARIABLES 52 

#*»**##**»»****** 53 

56 

NW - NUMBER OF POINTS WITH ZERO WEIGHT 55 

NS S - NUMBER OF SEGMENTS FOR NEW FIT (RETURNED FROM 56 

SUBROUTINE REFIT) 57 

DEV - DEVIATION OF FITTED CURVE FROM ORIGINAL DATA POINTS 58 

ERR - RELATIVE ERROR 59 
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•••••••••>•• 60 

SUMX - SUM OF Y « 61 

SUMX2 - SUM OF Y SQUARED • 62 

SUMY - SUM OF YC ** - USED IN CALCULATING 63 

SUMY2 - SUM OF YC SQUARED * CORRELATION COEFFICIENT 64 

SUMX Y - SUM OF Y TIMES YC » 65 

•••••*•••>•* 66 

************ 67 

SUMD - SUM OF DEVIATIONS *• - USED IN CALCULATING 68 

SUMD2 - SUM OF DEVIATIONS SQUARED * VARIANCE AND STANDARD 69 

•••»•*••**»• DEVIATION 70 

FN - NUMBER OF DEGREES OF FREEDOM 71 

VAR - VARIANCE 72 

STDEV - STANDARD DEVIATION 73 

CORR - CORRELATION COEFFICIENT 74 

CORMAX - MAXIMUM POSSIBLE CORRELATION COEFFICIENT 75 

76 

IM - MAXIMUM ORDER OF POLYNOMIALS <IM=3) 77 

78 

79 

80 

TO CHANGE THE MAXIMUM NUMBER OF POINTS OR THE MAXIMUM NUMBER DF 81 

SEGMENTS THE PROGRAM WILL FIT, THE FOLLOWING TWO VARIABLES MUST BE 82 

CHANGED - 83 

84 

I NX - MAXIMUM NUMBER OF DATA POINTS (INX IS NOW SET AT 350) 85 

IXM - MAXIMUM NUMBER OF SEGMENTS (IXM IS NOW SET AT 10) 86 

87 

THE FOLLOWING DIMENSIONED VARIABLES MUST BE CHANGED ALSO - 88 

89 

X, Y, W, XX, YY, WW, YC, AND NBLANK MUST HAVE DIMENSION INX 90 

XM MUST HAVE DIMENSION (IXM+1) 91 

LLOW AND LHIGH MUST HAVE DIMENSION IXM 92 

THE REMAINING ARRAYS MUST HAVE DIMENSIONS THAT CORRESPOND TO THE 93 

NUMBER OF SEGMENTS AND THE HIGHEST ORDER POLYNOMIAL - 94 

A(IXM,IM+1) XWX! IXM,IM+1,IM+1) 95 

YWXt IXM.IM+l) C ( IXM-1, IM, IM+1 ) 96 

B( IXM-l.IM, IM.IM) BB( IXM-1 , IXM-1 , IM, IM) 97 

98 

***,#»****»#***#***************#****#*-»**»***************#***»*#** 99 

100 

THE SUBROUTINE DUB I D IS NECESSARY FOR DOUBLE PRECISION OUTPUT ON 101 

THE LEWIS COMPUTER 102 

103 

104 

105 

DIMENSION TITLE! 12) ,FMT ( 12 ) . FMTM ( 12 > 106 

DIMENSION X! 350 ) , Y! 350 ) , W ( 350 ) , XX ( 3 50 ! , YY ( 350 ) , WW ( 3 50 ) , YC ( 3 50 ) , 107 

1 NBLANK ( 350 ) 108 

DIMENSION XM( 11 ) ,LLOW ( 10 ) , LHIGH! 10 1 109 

DIMENSION A(10,4),XWX(10,4,4),YWX!10,4),C(9,3,4),B!9,3,3,3), 110 

1 BB(9,9,3,3> 111 

DOUBLE PRECISION B, BB , A , YC , XWX , YWX , C 112 

DOUBLE PRECISION OE V, ER R , VAR , STDEV , CORR , CORMAX, SUMX , SUMY, SUMXY , 113 

1 SUMX2 , SUM Y 2 , SUMD , SUMD 2 114 

EXTERNAL DUBIO 115 

LOGICAL LREFI T, TRANX, TRANY.NPUNCH 116 

INX = 350 117 

IXM = 10 118 

IM = 3 119 

120 

SET DIMENSIONS OF ARRAYS IN SUBROUTINES 121 

122 

CALL ORD! X,Y,W, XX, YY.WW .NBLANK, 350 ) 123 

CALL SEG1XX.YY.XM, LLOW, LHIGH, 350, 11, 10) 124 

CALL DEF!XX,YY,WW,XM, LLOW, LHIGH, XWX, YWX, C , 350 , 1 1 , 10 , 4, 3 , 9 ) 125 

CALL SLV(C,XWX, YWX, A, B, BB»9»3,4,10 ) 126 

CALL RFTIXX, A, XM, LLOW, LHIGH, 350, 10, 11,4) 127 

CALL TRN(X,Y,XM, YC, 350, 11) 128 

CALL DXMIXX.XM, LLOW, LHIGH, 350, 11, 10) 129 

CALL DSPL1XX,XM,LL0W,LH IGH, 350, 11, 10) 130 

CALL DNSIXX.XM, LLOW, LHIGH, 350, 11, 10 ) 131 

CALL DHIXX.YY.XM, LLOW, LHIGH, 350, 11, 10) 132 

CALL OL(XX,YY,XM, LLOW, LHIGH, 350, 11, 10) 133 
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C 134 

KASES = 0 135 

80 WRITE (6,24) 136 

C 137 

C IF KASES = 0, READ A NEW TITLE AND A NEW SET OF X.Y.W DATA 138 

C IN ANY CASE, READ NEW VALUES FOR M , NX, NS, NP, MB , NF , REF IT, 139 

C TRANX, TRAN Y , TOL 140 

C 141 

IF {KASES. EQ.O) READ (5,4) TITLE 142 

81 READ (5,1) M,NX,NS,NB,NF»LREFIT»TRANX,TRANY»NPUNCH« TOL 143 

KASES = KASES— 1 144 

IF (KASES. GE.O) GO TO 84 145 

READ (5,4) FMT 146 

READ (5, FMT) { X ( I ) , Y( 1 1 , W ( I ) , I = 1 , NX ) 147 

C 148 

C IF ALL WEIGHTS ARE READ AS 0, SET ALL WEIGHTS TO 1 149 

C 150 

NW = 0 151 

DO 82 1=1, NX 152 

IF (W(I).LE.O.) NW = NW+1 153 

82 CONTINUE 154 

IF (NW.NE.NX) GO TO 84 155 

DO 83 1=1, NX 156 

83 W { I ) = 1. 157 

C 158 

C TEST INPUT TRIGGERS TO SEE IF ADDITIONAL DATA IS NEEDED TO DIVIDE 159 

C THE CURVE INTO SEGMENTS 160 

C 161 

84 IF (NB.NE.O) READ (5,4) FMTM 162 

IF (NB.NE.O) READ (5, FMTM) ( XM( I),I=1,NB) 163 

IF (N8.GT.IXM) NB= IXM 164 

IF (KASES.LT. 01 READ (5,1) KASES 165 

IF(TRANX.OR.TRANY) CALL TR ANSF ( X, Y , NX , XM, NB , TRANX , TR ANY ) 166 

WRITE (6,10) TITLE 167 

C 168 

C CHECK FOR SUFFICIENT NUMBER OF DATA POINTS AND CORRECT ORDER DF 169 

C POLYNOMIALS 170 

C 171 

IF (NX.GT.M.ANO. (M.EQ.2.0R.M.EQ.3) ) GO TO 85 172 

WRITE (6.23) M, NX 173 

GO TO 80 174 

C 175 

85 CALL OROER{ X,Y, W , NX ,XX, YY,WW, NXX ) 176 

CALL SEGMNT(XX, Y Y , X M, LL OW , LHI GH ,NXX , NS , NB , NF , M, TOL , IXM) 177 

CALL DEFMAT(XX, Y Y ,W W, XM , LLOW, LH IGH , NXX , NS , M, XWX , YWX.C) 178 

XMIIXM+1) = XX ( 1 ) 179 

CALL ASOLVE(C,XWX,YWX,A,NS,M) 180 

IF (NS.EQ.l) LREFIT=. FALSE. 181 

IF (LREFiri CALL REF I T ( XX, A , XM, LLOW, LHI GH, NXX ,NS, NSS , M, TOL ) 182 

C 183 

C WRITE OUTPUT DATA 184 

C WRITE ORDER OF POLYNOMIALS AND NUMBER OF SEGMENTS 185 

C 186 

86 WRITE (6,11) M, NS 187 

IF (NS. EQ.O) GO TO 80 188 

IF (M.EQ.3) GO TO 90 189 

C 190 

C WRITE EQUATION FOR FITTED CURVE FOR M=2 191 

C WRITE COEFFICIENTS IN STYLE FOR M=2 192 

C 193 

IF (TRANX. AND. TRANY) WRITE (6,39) 194 

IF (TRANX. ANO.. NOT. TRANY) WRITE (6,40) 195 

IF ( .NOT. TRANX. AND. TRAN Y) WRITE (6,41) 196 

IF ( .NOT. TRANX. AND. .NOT. TRANY) WRITE (6,42) 197 

WRITE (6,121 198 

WRITE (6,44) 199 

WRITE (6,13) (( A(N, J) , J = 1 , 3 ) . N=1 , NS ) 200 

GO TO 91 201 

C 202 

C WRITE EQUATION FOR FITTED CURVE FOR M=3 203 

C WRITE COEFFICIENTS IN STYLE FOR M=3 204 
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C 205 

90 IF (TRANX. AND. TRANY) WRITE (6,46) 2064 

IF (TRANX.AND..N0T. TRANY) WRITE (6,47) 207 

IF (. NOT. TRANX. AND. TRANY) WRITE (6,48) 208 

IF ( .NOT. TRANX. AND. .NOT. TRANY) WRITE (6,49) 209 

WRITE (6,12) 210 

WRITE (6,45) 211 

WRITE (6,14) (( A(N, J) , J = 1 , 4 ) , N=l, NS ) 212 

213 

CALCULATE FITTED VALUES OF Y AND BACK TRANSFORM DATA 214 

215 

91 DO 103 1=1, NX 216 

DO 100 N= 1 , NS 217 

NN = N 218 

IF (X(I).LE.XM(N)) GO TO 101 219 

100 CONTINUE 220 

101 VC(I) = A ( NN, M+ 1 ) 221 

DO 102 J=1.M 222 

JJ = M+l-J 223 

YC ( I ) = YC( I )*X( I )+A(NN, JJ) 224 

102 CONTINUE 225 

103 CONTINUE 226 

IF ( TRANX.OR.TR ANY) CALL 8 TRANS ( X , Y , XM, YC , NX, NS . TRANX, TRAN Y > 227 

228 

WRITE SPLINE JOINTS 229 

230 

WRITE (6,15) 231 

WRITE (6,16) XM( IXM+1 ), (XM( I ) , 1=1, NS) 232 

233 

234 

CALCULATE DEVIATION ANO RELATIVE ERROR 235 

CALCULATE SUMS FOR VARIANCE AND CORRELATION COEFFICIENT 236 

WRITE X,Y,Y*. DEVIATION, AND RELATIVE ERROR 237 

238 

SUMX = 0.000 239 

SUMY = O.ODO 240 

SUMXY = O.ODO 241 

SUMX2 = O.ODO 242 

SUMY2 = O.ODO 243 

SUMD = O.ODO 244 

SUM02 = O.ODO 245 

WRITE (6,21) 246 

DO 110 1=1, NX 247 

DEV = YC ( I I -Y ( I ) 248 

IF (Y( I I.NE.O.O) GO TO 111 249 

ERR = UEV/YCI I ) 250 

GO TO 112 251 

111 ERR = DEV/Y ( I I 252 

112 SUMX = SUMX+ Y ( I ) 253 

SUMY = SUMY+YC( I ) 254 

SUMXY = SUMXY+Y ( I ) • YC ( I ) 255 

SUMX2 = SUMX2+Y ( I )»Y{ I ) 256 

SUMY2 = SUM Y2+YC ( I) »YC( I) 257 

SUMD = SUMD+DEV 258 

SUMD2= SUM02+DEV*DEV 259 

11J WRITE (6,20) X( I),Y( I ) »YC{ I) , DEV, ERR 260 

261 

CALCULATE AND WRITE VARIANCE, STANDARD DEVIATION, AND CORRELATION 262 

COEFFICIENT 263 

264 

FN = FLOAT ( NX - M* ( NS -l ) ) 265 

FX = FLOAT ( NX ) 266 

VAR = (SUMD2-SUMD»SUMD/FX)/FN 267 

STDEV = SORT ( VAR ) 268 

CORR = FN»(FX»SUMXY-SUMX*SUMY)/FX/SORT( (FX*SUMX2-SUMX*SUMX)* 269 

1 ( FX*SUMY2— SUM Y*SUMY ) ) 270 

CORMAX = FN/FX 271 

WRITE (6,22) VAR, CORR, STDEV, CORMAX 272 

C 273 

IF (LREFIT) GO TO 88 274 

WRITE (6,43) 275 

GO TO 79 276 

88 IF (TOL.LT.O.) GO TO 78 277 

WRITE (6,37) 278 
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GO TO 79 279 

78 NS = NSS 280 

NB = NS— 1 281 

DO 77 1=1, NS 282 

NPTS = LHIGHl II 283 

XMU) = XX ( NPTS I 284 

77 CONTINUE 285 

IF (TRANX.OR.TR ANY) CALL TRANSF ( X . Y , NX , XM ,NB , TR AN X, TR ANY) 286 

CALL DEFMATIXX, Y Y, WW, XM ,LLOW, LHI GH , NXX , NS , M, XWX,YWX,C) 287 

CALL ASOLVE(C,XWX,YWX,A,NS,M) 288 

WRITE (6,38) 289 

WRITE (6,36) 290 

WRITE (6,24) 291 

WRITE (6,10) TITLE 292 

WRITE (6,36) 293 

L REF I T = .FALSE. 294 

GO TO 86 295 

79 IF (NPUNCH) GO TO 80 296 

WRITE (6,32) TITLE 297 

WRITE (6,33) NS, { XM (N ) , N= 1 , NS ) 298 

IF (M.EQ.3) GO TO 89 299 

WRITE (6,30) ( ( A (N, I) , 1=1,3) ,N=1,NS ) 300 

GO TO 80 301 

89 WRITE (6,31) ( ( A ( N. I ) , I = 1 , 4 ) , N= 1 , NS ) 302 

GO TO 80 303 

C 304 

1 FORMAT <5I3,4L3,F12.6) 305 

3 FORMAT (2413) 306 

4 FORMAT (12A6) 307 

10 FORMAT ( 1H ,5X,12A6) 308 

11 FORMAT ( 1H0 , 22HOEGREE OF POLYNOMIAL = , 1 5 , 10X , 20HNUMBER OF SEGMENTS 309 

1 = , I 5 ) 310 

12 FORMAT (1H0.41HSEGMENT COEFFICIENTS IN ASCENDING ORDER - ) 311 

13 FORMAT ( 1H ,1P3D25.15> 312 

14 FORMAT (1H .1P4D25.15) 313 

15 FORMAT (1H0.19HSPLINE JOINTS ARE -) 314 

16 FORMAT (IH ,21X,7G15.7> 315 

20 FORMAT ( 1H , 1P2 E 17 . 7, IP 2D25. 1 5, 0PD2 5. 1 5 ) 316 

21 FORMAT (1H0,9X, 1HX, 14X, 1HY , 23X, 2HY* , 2 3X , 3HDE V, 20X ,5HR— ERR ) 317 

22 FORMAT ( 1 HO , / , 1 HO , 50HCORREL AT I ON OF FITTED DATA TO ORIGINAL DATA 318 

1 ,/, 1H0 ,16X, 10HVARI ANCE = , 1P025. 15 , 12 X , 19HC0RREL AT I ON IND 319 

2EX =,0PD25. 15,/, IH ,6X, 20HSTANDARD DEVIATION = , 1PD2 5. 1 5, 1 OX , 320 

321HMAXI MUM CORRELATION =,OPD25.15> 321 

23 FORMAT ( 1 HO, 30H CANNOT MAKE VALID FIT WITH M =,I3,9H AND NX =,I3) 322 

24 FORMAT ( 1H1 ) 323 

30 FORMAT ( 1 H$ , 3D20 . 13 ) 324 

31 FORMAT ( 1H$,4D20.13) 325 

32 FORMAT (1HS.12A6) 326 

33 FORMAT ( IH* , I 3, ( / , 1 H$ , 5 E 14.7 ) ) 327 

36 FORMAT ( IH , 82HDUPL IC AT ION OCCURED IN FIRST SET OF COEFFICIENTS - 328 

1CURVE WAS REFIT IN NEW SEGMENTS) 329 

37 FORMAT ( IH ,43HN0 DUPLICATION IN FIRST SET OF COEFFICIENTS) 330 

38 FORMAT ( 1H0 , 20HR EF I T CHECK WAS MADE) 331 

39 FORMAT ( 1H0 , 69H E QUA T I ON FITTED IS LOG Y = AO + A1 (LOG X 332 

1) + A2 (LOG X ) » » 2 ) 333 

40 FORMAT (1H0 , 64HE QUAT I ON FITTED IS Y = AO + A1 (LOG X) + 334 

1A2 (LOG X ) **2 ) 335 

41 FORMAT ( 1HO.56HE0UATION FITTED IS LOG Y = AO + A1 X + A2 336 

1 X**2) 337 

42 FORMAT ( 1H0 , 52H EQUATI ON FITTED IS Y = AO + A1 X + A2 X»» 338 

12) 339 

43 FORMAT (1H0.19HN0 REFIT CHECK MADE) 340 

44 FORMAT ( 1H0 , 10X , 2HA0, 23X, 2HA1 , 23X, 2HA2 ) 341 

45 FORMAT ( 1 HO , 10X , 2HA0, 23 X, 2HA 1 , 23X , 2HA2 , 23X, 2H A3 ) 342 

46 FORMAT ( 1H0 , 84HEQUAT I ON FITTED IS LOG Y = AO ♦ A1 (LOG X 343 

1) "f A2 (LOG X ) • »2 + A3 (LOG X>**3) 344 

47 FORMAT ( 1H0 , 8 IH E QUA T I ON FITTED IS Y = AO + A1 (LOG X) + 345 

1A 2 (LOG X ) **2 ♦ A3 (LOG X)**3 ) 346 

48 FORMAT (1 HO, 66H EQUATI ON FITTED IS LOG Y = AO + A1 X + A2 347 

1 X**2 + A3 X»»3 ) 348 

49 FORMAT ( 1H0 , 63H EQUA TI ON FITTED IS Y = AO + A1 X «■ A2 X** 349 

12 ♦ A3 X*»3 ) 350 

C 351 

END 352 
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APPENDIX D 


VARIABLES USED BY SEVERAL SUBROUTINES 

The variables used by several subroutines of the program FITLOS are defined as 
follows: 

XM Array of spline joints. 

LLOW Array of indices of the first point in each subset. LLOW(l) = 1 and 

LLOW(N) = lowest I such that XM(N - 1) < X(I) < XM(N) for 
N = 2, .... NS. 

LHIGH Array of indices of the last point in each subset. LHIGH(N) = highest I 

such that XM(N) < X(I) < XM(N + 1) for N = 1, . . . , NS - 1, and 
LHIGH (NS) = NX. 

T -1 

XWX Multidimensioned array (X WX) . The subscripts on XWX have the 

T “1 

same order as the subscripts on matrix (X WX) of appendix B. 

T 

YWX Multidimensioned array Y WX. The subscripts correspond to the sub- 

T 

scripts on vector Y WX of appendix B. 

C Multidimensioned array of constraints. The subscripts correspond to the 

subscripts on matrix C of appendix B. 

A Multidimensioned array of undetermined coefficients. The subscripts 

correspond to the subscripts of vector A of appendix B. 

X Array of the ordered independent variable. 

Y Array of dependent variable that corresponds to X. 

In the main program FITLOS, X and Y are the names of the input arrays, while 
XX and YY are the names of the ordered data. 
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APPENDIX E 


DESCRIPTION OF SUBROUTINES 

The subroutines of the program FITLOS are described in this appendix. After the 
descriptions, all the subroutines are listed followed by all the flow charts. 

TRANSF 

Subroutine TRANSF makes a base 10 log transformation on the input data. If any of 

-30 

these data are not greater than zero, TRANSF changes that number to 10 
BTRANS 

Subroutine BTRANS converts the transformed data back to its original form. Since 
YC, the calculated values of y*, are actually log^YC), these data are also back trans- 
formed so they have the same form as the input data. 

ORDER 

Subroutine ORDER arranges the input data in order of ascending x. Since the defi- 
nition of a spline function requires that x^ < x i+ ^, ORDER averages the y's for which 
duplicate x's occur. This average y is a weighted average, 


y = 


E y.w. 

J 

j 


2>j 


The total weight, X/ W j> 


becomes the weight of the average point. 


j 

data, the ordered data are put into new arrays. 


To preserve the input 


SEGMNT 

Subroutine SEGMNT determines the spline joints and the low and high indices of the 
points in each subset. SEGMNT first tests the variable NB. If NB # 0, the spline 
joints have been supplied by the user. In that case, SEGMNT calls subroutine DIVXM to 
determine the index arrays. If NB = 0, SEGMNT then tests the variable NS. 

If NS + 0, the number of segments has been chosen by the user. In that case, 
SEGMNT calls subroutine DIVNS to divide the data as evenly as possible among the NS 
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subsets. Subroutine DIVNS calculates the index arrays and determines the spline joints. 

If NS = 0, SEGMNT tests the number of data points NX. 

If NX < 3M, SEGMNT calls subroutine SPESHL to make a special division of the 
data into one or two subsets with specially determined spline joints. These special 
spline joints and the index arrays are listed in the main-text section HOW DATA ARE 
DIVIDED INTO SUBSETS. 

If NX > 3M, SEGMNT sets the number of subsets to be the maximum possible num- 
ber based on the number of data points and the degree of the polynomial. This number 
is (NX - 1)/M. 

SEGMNT then tests the variable NF. If NF < 0, subroutine FFLOW is called to 
do a force -fit division starting at the low end of the data. If NF > 0, subroutine FFHIGH 
is called to do a force-fit division starting at the high end of the data. If NF = 0, sub- 
routine DIVNS is called with NS = (NX - 1)/M. 

DIVXM 

Subroutine DIVXM divides the data into subsets according to spline joints (xm) chosen 
by the user. DIVXM first puts the (xm) in ascending order. Then it eliminates any of 
the (xm) that are outside the range of x and adjusts the number of spline joints NB 
accordingly. 

DIVXM then determines the indices of the first and last points in each subset. Then 
it checks whether each subset has a sufficient number of points. If LHIGH(I) - LLOW(I) 

+ 1 < M, there are not enough points in subset I and that subset must be combined with 
its neighbors. DIVXM also changes the spline joints to correspond to the new index 
arrays. 

DIVNS 

Subroutine DIVNS divides the data into NS subsets as evenly as possible. DIVNS 
first makes NS a "proper” number. It chooses the smallest of three possible values 
which are as follows: the chosen NS, the maximum number of subsets based on the 
number of data points and the degree of the polynomial, and the dimension of the arrays 
LLOW and LHIGH, which is called LEM in the program. 

In dividing the data as evenly as possible, DIVNS uses fixed-point arithmetic to 
eliminate the possibility of a fractional number of points in a subset. The spline joints 
and the index arrays are determined as the division takes place. 

SPESHL 

Subroutine SPESHL makes an arbitrary division of the data into one or two subsets. 

If the number of points NX is not greater than 2M, only one segment is possible. If 
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NX is between 2M and 3M, SPESHL divides the data into two subsets where the spline 
joints and index arrays are defined in the main -text section HOW DATA ARE DIVIDED 
INTO SUBSETS. 

FFLOW 

Subroutine FFLOW divides the data into subsets by force -fitting starting at the low 
end of the data. The maximum possible number of subsets NS appears in the calling 
vector. FFLOW first sets the index arrays LLOW and LHIGH to zero. It then starts 
force -fitting as described in the main body of the report. For a point to be accepted in 
a subset it must fall on the Lagrange polynomial within a given amount of precision 
TOL; that is, the ratio | y (cal c)/y (given) j = 1 + TOL. For y(given) = 0, the acceptance 
criterion is |y(calc) - y(given) | < TOL. Spline joints are determined as the last point 
in each subset. 

FFHIGH 

Subroutine FFHIGH does a force -fit division of the data into subsets starting at the 
high end of the data. It first sets the index arrays to zero. Then it starts force -fitting, 
but begins with the NS^ 1 segment. Consequently, some low -order elements of LLOW 
and LHIGH could remain zero. If they do, the nonzero elements are moved down so 
that LLOW(l) = 1. The elements of LHIGH and XM are moved down simultaneously 
and the value of NS is reduced accordingly. In all other respects, however, FFHIGH 
and FFLOW are essentially the same. 

REFIT 

Subroutine REFIT checks whether the curve was fit in more segments than were 
necessary. To do this, it checks whether the coefficients for a low -order segment would 
give the same value of y* for points in a higher order subset as the coefficients for the 
higher order subset. Subroutine REFIT works in essentially the same way as the force- 
fitting subroutines except the test polynomial is defined by the coefficients from the lower 
order segment instead of a Lagrange polynomial. The use of TOL is the same as in 
subroutine FFLOW. 

MINVRT 

Subroutine MINVRT inverts a double -precision matrix by Gaussian elimination 
(ref. 6). It also calculates the determinant of the matrix. If the determinant is zero, 
that is, if the matrix is singular, an error message is printed and the null matrix is 
returned to the calling program. If the matrix is nonsingular, MINVRT finishes the 
Gaussian elimination. Pivoting is not necessary since the matrices are small and well 
conditioned. Then the inverse is multiplied by the input matrix and the maximum devia- 
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tion of the elements of the product matrix from the elements of the identity matrix is 
returned to the calling program. This measures the accuracy of the inverse. Finally, 
the inverse is transferred to the input matrix for return to the calling program. 

DEFMAT 

Subroutine DEFMAT defines the matrices (X T WX) -1 , Y T WX, and C from the 
arrays of ordered data and the array of spline joints. The multiple subscripts on the 
arrays XWX, YWX, and C correspond to the subscripts on matrices (X^WX)"\ 

Y^WX, and C of appendix B. 

ASOLVE 

Subroutine ASOLVE solves equation (5). If there is only one segment, the simple 
matrix multiplication A = (X T WX)" 1 (Y T WX) 1 is performed. For more than one seg- 
ment, matrix B is defined by equation (6) of appendix B. Since each row of B has only 
three nonzero submatrices, only these three submatrices are calculated. Then B is 
inverted by the process described in appendix B. The E and DELTA matrices are 
the same as the E and A matrices defined in appendix B. 

Beginning with statement 500, the matrix multiplication of equation (5) is performed. 
Since there are four types of elements in the matrix product C T B -1 C, there are four 
separate techniques used for calculating these elements. These four types of elements 
are defined in appendix B . When the matrix product C T B -1 C has been formed, the re- 
maining multiplication is finished. The vectors V and VV are the same as defined in 
appendix B. 


39 


oon nonono 


$ I BFTC TR1 

1 

PROGRAM VARIABLES 2 

A 

YC - CALCULATED VALUES OF Y 5 

6 

SUBROUTINE TRlVI X, Y, XM, YC, IX, IXM) 7 

DIMENSION X ( IX) , Y( I XI , YCI IXI , XM I IXM) 8 

DOUBLE PRECISION YC 9 

GO TO 3 10 

11 

TRANSFORMATION SUBROUTINE 12 

13 

ENTRY TRANSFIX, Y , NX , X M , NB , T RANX , TR ANY ) 14 

LOGICAL TRANX.TRANY 15 

IF I.NOT.TRANX) GO TO 1 16 

DO 10 1=1 , NX 17 

IF (X(I).LE.O.) XII )= 1, E-30 18 

10 XIII = ALOGIOIXI I ) ) 19 

IF (N8.E0.0) GC TO 1 20 

DO 11 1=1, NB 21 

IF (XM(I).LE.O.) XM ( I ) = 1 . E— 30 22 

11 XM ( I ) = ALOGIOI XM I I ) ) 23 

1 IF I . NOT. TR ANY ) GO TO 3 2A 

DO 12 1=1, NX 25 

IF (Y(I).LE.O.) Y ( I 1 = 1. E-30 26 

12 YII) = AL0G10IYI II) 27 

GO TO 3 28 

C 29 

ENTRY BTRANSIX, Y.XM.YC, NX, NS, TRANX.TRANY) 30 

IF I . NOT. TRANX) GC TO 2 31 

DO 13 1=1, NX 32 

13X11) = 10.»*X(I ) 33 

00 14 1=1, NS 34 

14 X M I I ) = 10. **XM I I I 35 

XM(ll) = 10 . »*XM 111) 36 

2 IF I . NOT . TRANY) GC TO 3 37 

DO 15 1=1, NX 38 

YII) = 10.**Y(I) 39 

16 YC I I) = 10 . * * YC I I ) 40 

C 41 

3 RETURN 42 


tlBFTC ORDR 

C 1 

C PROGRAM VARIABLES 2 

C #*«,,»*«*,,*»*#** 3 

C A 

C XT - ORIGINAL VALUES OF THE INDEPENDENT VARIABLE 5 

C YT - ORIGINAL VALUES OF THE DEPENDENT VARIABLE 6 

C WT - ORIGINAL WEIGHTS 7 

C NXT - NUMBER OF ORIGINAL POINTS 8 

C 9 

C X - ORDERED ARRAY OF INDEPENDENT VARIABLFS 10 

C Y - ORDERED ARRAY OF DEPENDENT VARIABLES 11 

C W - ORDERED ARRAY OF WEIGHTS 12 

C NX - NUMBER OF ORDERED DATA POINTS 13 

C 14 

C NBLANK - BOOKKEEPING ARRAY, N8LANK (11=0 MEANS POINT I HAS 15 

C BEEN TRANSFERED TO THE NEW ARRAYS 16 

C KK - INDEX OF THE AVERAGED POINT IN THE NEW ARRAYS 17 

C N - NUMBER CF POINTS WITH SAME XT VALUE 18 

C SUMY - SUM OF YT VALUES FOR POINTS WITH SAME XT VALUE 19 

C SUMW - SUM OF WEIGHTS FOR POINTS WITH SAME XT VALUE 20 

C 21 

SUBROUTINE ORD(XT,YT,WT,X,Y,W, NBLANK, IX) 22 

DIMENSION XT< IX ) , YT ( I X) ,WT( IX) ,X( IX) , Y( IX) ,W( I X) .NBLANK! IX) 23 

GO TO 300 24 
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25 

ARRANGE DATA IN ORDER OF ASCENDING X AND AVERAGE Y FOR WHICH 26 

DUPLICATE VALUES OF X OCCUR 27 

28 

ENTRY ORDERIXT, YT , W T, NX T , X, Y , W , NX > 29 

30 

DEBUG ( XT ( I) , YT (I) , WT I I) , 1 = 1, NXT) 31 

DO IOC 1=1, NXT 32 

100 NBLANKl I ) = 1 33 

NX = 0 34 

KK = 1 35 

DO 230 I = 1 ,NX T 36 

IF ( NBLANKl I ) .EQ.O) GO TO 230 37 

N = 1 38 

SUMY = YT 1 I ) *WT (I) 39 

SUMW = WT ( I) 40 

NBLANKl I ) = 0 41 

11=1+1 42 

DO 200 J= I I * NXT 43 

IF (NBLANKl Jl.EQ. 01 GO TO 200 44 

IF I XT ( J) • NE . XT I I ) ) GO TO 200 45 

SUMY = SUMY+YTI J ) »WT1 J ) 46 

SUMW = SUMW+WTl J ) 47 

N = N+l 48 

NBLANK(J) =0 49 

200 CONTINUE 50 

IF (KK.EO.l) GU TO 221 51 

00 220 J= 1 , NX 52 

IF (X(J).LE.XT(I)) GO TO 220 53 

KN = NX — J + 1 54 

DO 210 K=1,KN 55 

KK = NX+2-K 56 

X ( KK ) = X(KK-l) 57 

Y ( KK I = Y(KK-l) 58 

21' W ( KK ) = w(KK-l) 59 

GO TO 222 60 

22u CONTINUE 61 

222 KK = KK— 1 62 

221 X(KK) = X T ( I ) 63 

Y ( KK ) = SUMY/ SUMW 64 

W ( KK ) = SUMW 65 

NX = NX+1 66 

KK = NX +2 67 

23C CONTINUE 68 

DEBUG ( X< I ) . Y ( I ) ,W( II , 1=1, NX) 69 

70 

30' RETURN 71 

END 72 


BFTC SGMNT 

SUBROUTINE S EG( X , Y , XM , L LOW , LH I GH, I X , I XM , I L ) 1 

DIMENSION X( IX) , Y( IX) ,XM( IXMJ.LLOWI IL).LHIGH(IL) 2 

RETURN 3 

4 

OIVIDE DATA INTO SUBSETS BY DETERMINING SPLINE JOINTS AND 5 

THE NUMBER OF POINTS IN EACH SUBSET 6 

7 

ENTRY SEGMNT (X, Y, XM.LLOW, LHIGH, NX, N S , NR , NF , M , TOL ,LI M) 8 

9 

DIVIDE ACCORDING TO PREDETERMINED BREAK POINTS 10 

11 

200 IF INB.EO.O) GO TO 400 12 

CALL DIVXM(X,XM,LLOW,LHIGH,NX,NS,NB,M) 13 

WRITE (6,21) 14 

RETURN 15 

16 

DIVIDE ACCORDING TO PREDETERMINED NUMBER OF SEGMENTS 17 

18 

400 IF (NS.EO.O) GO TO 500 19 

CALL DIVNStX, XM.LLOW, LHIGH, NX, NS, NB.M.LIM) . 20 
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WRITE (6,23) 21 

RETURN 22 

C 23 

C DO THE NUMBER OF POINTS REQUIRE A SPECIAL DIVISION 24 

C 25 

500 IF INX.GT.3«M) GO TO 600 26 

CALL SPESHL!X,XM,LLOW,LHIGH,NX,NS,NB,M) 27 

WRITE (6,27) NX 28 

RETURN 29 

C 30 

C DIVIDE ACCORDING TO FORCE FIT SCHEME OR AS EVENLY AS POSSIBLE 31 

C AMONG SEGMENTS 32 

C 33 

600 NS = (NX-1)/M 34 

IF INF) 610,620,630 35 

610 CALL FFLOW( X,Y, XM,LLOW, LHIGH.NX.NS, NB, M, TOL ) 36 

WRITE (6,24) 37 

RETURN 38 

C 39 

620 CALL DIVNS(X,XM,LLOW,LHIGH,NX,NS,NB,M,LIM) 40 

WRITE (6,25) 41 

RETURN 42 

C 43 

630 CALL FFHIGH(X,Y ,XM,LLOW,LHI GH ,NX,NS,NB»M,TOL) 44 

WRITE (6,26) 45 

RETURN 46 

C 47 

20 FORMAT ( 1 HO .47HNUMB ER OF DATA POINTS REQUIRES SPECIAL DIVISION) 48 

21 FORMAT ( 1H0.34HSPLINE JOINTS CHOSEN BY PROGRAMMER) 49 

22 FORMAT ( 1H0 , 52HNUMB ER OF POINTS IN EACH SUBSET CHOSEN BY PROGRAMME 50 

1R) 51 

23 FORMAT (1H0, 83HDATA DIVIDED AS EVENLY AS POSSIBLE AMONG THE NIUMBE 52 

1R OF SUBSETS CHOSEN BY PROGRAMMER) 53 

24 FORMAT (1H0.74HDATA DIVIDED INTO SUBSETS BY FORCE FITTING STARTING 54 

1 AT THE LOW END OF DATA) 55 

25 FORMAT ( 1H0 , 70HDATA DIVIDED AS EVENLY AS POSSIBLE AMONG THE MAXIMU 56 

1M NUMBER OF SUBSETS) 57 

26 FORMAT ( 1H0 , 75HDATA DIVIDED INTO SUBSETS BY FORCE FITTING STARTING 58 

1 AT THE HIGH END OF DATA) 59 

27 FORMAT ( 1 HO , 1 5, 2X , 3 2HP0 IN TS REQUIRES SPECIAL OIVISION) 60 

C 61 

END 62 


IBFTC DVXM 

1 

PROGRAM VARIABLES 2 

3 

4 

T - TEMPORARY STORAGE USED IN ORDERING SPLINE JOINTS 5 

KST - INDEX OF FIRST POINT IN NEW SUBSET 6 

NSS - SUBSET COUNTER WHEN A SUBSET DOES NOT HAVE SUFFICIENT 7 

POINTS 8 

NPLUS - NUMBER OF POINTS IN (1*1) SUBSET 9 

NP2 - ONE HALF THE NUMBER OF POINTS IN SUBSET I 10 

11 

SUBROUTINE DXM( X , XM ,L LOW, LH I GH, I X, I XM , I L ) 12 

DIMENSION XI IX) ,XM( IXM) ,LLOW( IL ) ,LHIGH( IL ) 13 

RETURN 14 

15 

DIVIDE ACCORDING TO PRECHOSEN SPLINE JOINTS 16 

17 

ENTRY DIVXM(X,XM,LLOW,LHIGH,NX,NS,NB,M) 18 

19 

20 

CHECK THAT SPLINE JOINTS MATCH THE RANGE OF X 21 

22 

DEBUG ( XM( I ) , 1=1 ,NB) 23 

300 DO 310 1=1, NB 24 

IF ( I .EO.NB) GO TO 310 25 
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II = 1 + 1 26 

DO 305 J= 1 1 , NB 27 

IF (XM(I).LE.XM(J)) GO TO 305 28 

T = XM( I ) 29 

XM ( I ) = XM(J) 30 

XM(J) = T 31 

305 CONTINUE 32 

310 CONTINUE 33 

DEBUG IXM(I J,I=1,NB> 36 

C 35 

311 DO 320 1=1, NB 36 

IF t I.E0.NB1 GO TO 320 37 

IF IXMI I).GE.XtH+l) ) GO TO 320 38 

II = 1+1 39 

DO 315 J= 1 1 , NB 40 

315 XM(J-l) = XM(J> 41 

NB = NB— 1 42 

SO TO 311 43 

320 CONTINUE 44 

DEBUG NB, CXMI I) , I=1,NB) 45 

C 46 

NS = NB+1 47 

DO 330 1=1, NB 48 

IF IXM( I ) .LT.XINX) ) GO TO 330 49 

NS = NS— 1 50 

330 CONTINUE 51 

NB = NS— 1 52 

XM(NS) = X ( NX ) 53 

DEBUG NS, <XM< I) , 1=1 ,NSI 54 

55 

DETERMINE LOW AND HIGH INDICES 56 

57 

L LOW (11= 1 58 

X ST= 1 59 

DO 350 1=1, NS 60 

DO 340 K=KST,NX 61 

IF (XM(I).GT.X(KI) GO TO 340 62 

LHIGHI I ) = K-l 63 

IF(K.EO.NX) GO TO 340 64 

IF IXMIII.EO.XIKn LH I GUI I ) =K 65 

L LOW ( 1+1) = LHIGHI I ) 66 

IF IXMf I).NE.X(K) ) LLOWI I+1)=LHIGH( I ) + 1 67 

X S T= K + 1 68 

SO TO 350 69 

340 CONTINUE 70 

350 CONTINUE 71 

LHIGHINS) = NX 72 

DEBUG ( L LOW (II, 1=1, NS) 73 

DEBUG ( LHIGHI I ), 1=1 , NS) 74 

75 

CHECK FOR SUFFICIENT POINTS IN EACH SUBSET 76 

77 

11=0 78 

DO 360 1=1, NS 79 

I FILHIGHI I l-LLOW I I ) +1.GT.M) GO TO 360 80 

11= I 81 

WRITE (6, 10) I 82 

10 FORMAT ( 1H0.29H INSUFF IC IENT POINTS IN SUBSET, 15) 83 

360 CONTINUE 84 

85 

IF ANY SUBSETS ARE DEFICIENT, COMBINE THEM WITH OTHER SUBSETS 86 

87 

IF(II.EQ.O) RETURN 88 

N SS= NS 89 

400 DO 470 1=1, NSS 90 

NPTS= LHIGHI I >-LLOW( I )+ 1 91 

IFINPTS.GT.M) GO TO 470 92 

IFII.NE.l) GO TO 410 93 

L H I GH I 1 ) = LHIGHI2) 94 

L H IGH I 2 ) =0 95 

L LOW I 2 ) =0 96 

XM ( 1 ) = XM ( 2 ) 97 

DEBUG I 98 

GO TO 480 99 
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410 IFII.NE.NSS) GO TO 420 100 

LHlGHINSS-l) = NX 101 

L HIGH ( NSS ) =0 102 

LLOW { NS S ) =0 103 

XM ( NSS- 1 ) = X(NX) 104 

DEBUG 1,1 105 

GO TO 480 106 

420 NPLUS= LHIGHI 1+1 l-LLOWI 1+1 )+l 107 

IF(NPLUS.GT.M) GO TO 430 108 

LHIGHI I )= LHIGH1 1+1 ) 109 

L HIGH 1 1 ♦ 1 ) = 0 110 

LLOW( I + 1 ) =0 111 

XHI1 ) = XM( 1+1) 112 

DEBUG 1,1,1 113 

SO TO 480 114 

430 NP2 = (LHIGHI D-LLOWI I)+l)/2 115 

IF (NP2»2. NE. LHIGHI I)-LLOWt I )+l) GO TO 460 116 

LHIGHI I -1 ) = LHIGHI 1-1 )+NP2 117 

L HIGH ( I ) =0 118 

L LOW ( I + 1 ) = LLOW I 1 + 1 1-NP2 119 

IF (LLOWI I + D.LT. LHIGHI I-m LLOW 1 1 + 1 ) =LH I GH I I- 1 ) 120 

L LOW I I ) =0 121 

DEBUG I , I , I , I 122 

440 I FI LLOWI 1+1 ) .NE. LHIGHI 1-1 )) GO TO 450 123 

NP2 = LLOWI I + 1 ) 124 

X M I 1-11= XINP2) 125 

DEBUG I , I , I. I, I 126 

GO TO 480 127 

450 NP2 = LLOWI 1*1 ) 128 

XMII — 11= XINP2) 129 

NP2- LHIGHII-l) 130 

XM (1-11= . 5« I XM I 1-1 ) + X I NP 2 ) ) 131 

DEBUG I, I , I, I, I, I 132 

GO TO 480 133 

460 LHIGHI I— 1 )= LHIGH11-11+NP2 134 

LHIGH I I ) =0 135 

LLOWI I + 1 ) = LLOWI I+11-NP2 136 

LLOW I I ) -0 137 

GO TO 440 138 

470 CONTINUE 139 

NS = NSS 140 

RETURN 141 

142 

COMPACT INDEX ANO SPLINE JOINT ARRAYS AND CHECK AGAIN 143 

144 

480 00 500 1=1, NSS 145 

DEBUG I , LLOWI I ), LHIGHI I 1 146 

IFILLnwl I ) .GT.O. AND. LHIGHI I ) .GT.O) GO T0500 147 

11= I 148 

NST = NSS-l 149 

DO 490 J= I I , NST 150 

LLOWI J)= LLOWIJ+l) 151 

LHIGHI J )= LHIGHI J+l I 152 

XMI J)= XMI J+l) 153 

490 CONTINUE 154 

NSS = NSS-l 155 

DEBUG NSS 156 

DEBUG I LLOWIJ ), J=1,NSS) 157 

DEBUG 1 LHIGHI J) , J=1 , NSS I 158 

DEBUG I XMI J ) , J= 1 ,NSS) 159 

GO TO 400 160 

500 CONTINUE 161 

NS = NSS 162 

DEBUG I XMI I ), 1=1, NS ) 163 

DEBUG I LHIGHI I ), 1 = 1 , NS) 164 

DEBUG I LLOW (I) , 1=1, NS) 165 

C 166 

RETURN 167 

END 168 
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IBFTC DVNS 

PROGRAM VARIABLES 

»*«»*****»»•**»•» 

NSCRIT - MAXIMUM NUMBER OF SEGMENTS BASED ON DEGREE OF 
POLYNOMIAL AND NUMBER OF DATA POINTS 
NS - SMALLEST OF THE THREE POSSIBLE NUMBER OF SEGMENTS 
NPLFT - NUMBER OF POINTS THAT HAVE NOT BEEN ASSIGNED TO A 
SUBSET 

NSLFT - NUMBER OF AVAILABLE SUBSETS 
I - SUBSET INDEX 

NPTS - NUMBER OF POINTS THAT WILL BE IN THE I(TH) SUBSET 
LIM - DIMENSION OF ARRAYS LLOW AND LHIGH IN MAIN PROGRAM 

SUBROUTINE DNSI X , XM ,L LOW , LH I GH, IX, IXM, IL I 
DIMENSION X< IX) , XM< IXM) ,LLOW( I L ) , LH IGH ( I L ) 

GO TO 610 

DIVIDE ACCORDING TO PREDETERMINED NUMBER OF SEGMENTS 

ENTRY DIVNS(X,XM,LLOW,LHlGH.NX,NS,NB,M, ILL) 

IF (NS.EQ.O) GO TO 600 

ONE SEGMENT REQUIRES SPECIAL HANDLING 

IF (NS.NE.l) GO TO 500 
L LOW l 1 ) = 1 
L HIGH ( 1 ) = NX 
NB = 0 
RETURN 

MORE THAN ONE SEGMENT MEANS DIVIDING THE DATA AS EVENLY AS 
POSSIBLE AMONG THE SEGMENTS 

500 N SCR I T = (NX-1) /M 

510 NS = MI NO (NS.NSCRIT ,L IM ) 

DEBUG NS.NSCRIT, LIM 

511 NPLFT = NX 
NSLFT = NS 
NN = NS-1 
LLOW(l) = 1 
DO 520 1=1, NN 

NPTS = NPL FT/NS LFT+1 

NSLFT = NSLFT-1 

NPLFT = NPLFT-N P T S+ 1 

LHIGH(I) = LLOW I I ) +NPTS-1 

IF (I.LT.NS) LLOW ( I +1 I =LHIGH(I) 

NPTS =L H I GH ( I ) 

X M ( I ) = XINPTS) 

520 CONTINUE 

L H I GH (NS) = NX 
XM(NS) = X ( NX ) 

522 NB = NS-1 

DEBUG NS , NB 

DEBUG I LLOW ( I ) ,1 = 1, NS) 

DEBUG ( LH IGH ( I) , 1=1, NS) 

DEBUG (XM (I), 1=1, NS) 

610 RETURN 


600 WRITE (6,10) 

10 FORMAT ( 1 HO , 1 5HNS = 0 IN DIVNS) 
RETURN 
END 
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JIBFTC SPSHL 

SUBROUTINE DSPL < X , X M, LLOW , LHI GH, I X , IX M, I L ) l 

DIMENSION X( IX) ,XM{ IXM) ,LLOWt IL ) .LHIGHIIL ) 2 

RETURN 3 

C 4 

C SPECIAL DIVISION INTO SEGMENTS WHEN NUMBER OF POINTS IS BETWEEN 5 

C M+l AND 3M 6 

C 7 

ENTRY SPESHLIX, XM , L LOW , LH I GH, NX , NS . NB , M ) 8 

C 9 

C NUMBER OF POINTS LESS THAN 2M+1 REQUIRES ONE SEGMENT 10 

C 11 

100 IF <NX.GT.2«M) GO TO 200 12 

NS = 1 13 

NB = 0 14 

XM( 1 ) = X(NX) 15 

LLOW ( 1 ) =1 16 

LHIGH(l) = NX 17 

RETURN 18 

C 19 

200 NS = 2 20 

NB = 1 21 

GO TO ( 1,2.3) . M 22 

C 23 

1 WRITE (6,10) 24 

10 FORMAT (1H0.13HM CANNOT BE 1) 25 

CALL EXIT 26 

C 27 

C M = 2 28 

C 29 

2 LLOW ( 11=1 30 

L LOW ( 2 ) =3 31 

L H IGH ( 2 ) = NX 32 

XM (2 ) = X ( NX ) 33 

IF (NX.GT.5) GO TO 210 34 

XM ( 1 ) = X ( 3 ) 35 

LHIGH(l) =3 36 

RETURN 37 

C 38 

210 X M ( 1 ) = . 5*1 X(3 ) + X( 4) ) 39 

LHIGH(l) =4 40 

RETURN 41 

C 42 

C M = 3 43 

C 44 

3 L LOW ( 1 ) = 1 45 

X M ( 2 ) = X ( NX ) 46 

LHIGHI2) = NX 47 

IF (NX.GT.8) GO TO 310 48 

IF (NX.GT.7) GO TO 300 49 

X M ( 1 I = X<4) 50 

LLOW ( 2 ) =4 51 

LHIGH(l) =4 52 

RETURN 53 

C 54 

300 XM ( 1 ) = • 5* ( X ( 4 ) +X ( 5) ) 55 

L LOW ( 2 ) =5 56 

L HIGH ( 1 ) =4 57 

RETURN 58 

C 59 

310 LLOW ( 2 ) =5 60 

LHIGH(l) =5 61 

XM ( 1 ) = X ( 5 ) 62 

RETURN 63 

END 64 
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tIBFTC FFLW 

1 

PROGRAM VARIABLES 2 

3 
A 

NN - TRIAL NUMBER OF SEGMENTS 5 

NS - NEW SEGMENT COUNTER 6 

10 - INDEX OF FIRST POINT USED TO DETERMINE LAGRANGE 7 

POLYNOMIAL 8 

11 - INDEX OF SECOND POINT 9 

12 - INDEX OF THIRD POINT - LAST POINT FOR A QUADRATIC 10 

13 - INOEX OF LAST POINT FOR CUBIC 11 

NST - INDEX OF FIRST POINT TO BE TESTED 12 

* 13 

A » 1 A 

B *- INTERMEDIATE VALUES TO SIMPLIFY CODING OF THE 15 

C * LAGRANGE POLYNOMIAL 16 

0 * 17 

**••»** 18 

YJ - Y AT X(J> EVALUATED BY THE LAGRANGE POLYNOMIAL 19 

20 

SUBROUTINE DL ( X . Y , XM, LLOW , LH IGH, I X , I XM , I L I 21 

DIMENSION X< I X) , Yt I X) ,XM< IXM) ,LLOW( IL I , LHIGHI IL I 22 

RETURN 23 

2 A 

DETERMINE SUBSETS BY FORCE FITTING STARTING AT LOW END OF DATA 25 

26 

ENTRY FFLOW(X,Y,XM,LLOW.LHIGH,NX,NS,NB,M,TOL) 27 

28 

700 JJ = 1 29 

NN = NS 30 

NS = 0 31 

DO 710 1=1, NN 32 

LLOW ( I ) =0 33 

LHIGHI I I = 0 3 A 

710 CONTINUE 35 

LLOW 111= 1 36 

DO 750 N= 1 , NN 37 

10 = JJ 38 

NS = NS+1 39 

11 = 10 + 1 AO 

12 = 10 + 2 A 1 

13 = 10+3 A2 

NST = IO+M+1 43 

DEBUG NS, NST, 10, XI 10) AA 

IF (NST.GT.NX-M) GO TO 760 A5 

LHIGH(N) = NST- 1 A6 

DO 7A0 J=NST , NX A7 

JJ = J-l A8 

A = (XI Jl-XI 10) ) /(X 1 1 2 J —X I 1 1) ) A9 

B = I XI J)-X( 1 1) ) /IX (I2I-XI 10) ) 50 

C = (X(J)-X(I2) )/(X(Il)— X(IO) ) 51 

IF (M.EQ.3) GO TO 720 52 

YJ = Y( IO)*B»C-YI II )*C*A+Y( I2)*A»B 53 

DEBUG A,B,C,YJ,Y( J) ,X( J ) 5A 

GO TO 730 55 

720 D = XI J)-XI 13) 56 

YJ = D»(-Y( 1 0 ) » B*C/ IXI I 3 )— X I 10) )+Y! ll)»C*A/(X(I3)— XI ID 1 — Y l 12) » 57 

1 A*B/(X(I3)-XII2)))+Y(I3)»(XIJ)-XII0))*(XIJ)-XII1))*(X(J)-XII2))/ 58 

2 (XI I 3)-X( 10) )/ I XI I3)-X (I 1) )/ (XI I3I-XI 12) ) 59 

DEBUG A,B,C,D,X I J ) , Y( J ) , YJ 60 

730 IF (Y(J).EQ.O.) GO TO 731 61 

IF I ABS I l.-YJ/Y I J ) ) .GT. TOL ) GO TO 735 62 

LHIGHIN) = J 63 

GO TO 7 AO 6A 

731 IF I ABSIYJ-YI J) ) .GT.TOL ) GO TO 735 65 

L H I GH I N ) = J 66 

7A0 CONTINUE 67 

735 IF (N.NE.NN) L LOW (N+ 1 )=LH I GH I N ) 68 

750 CONTINUE 69 

760 L H I GH (NS) = NX 70 

71 

SELECT SPLINE JOINTS 72 
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C 73 

N8 = NS-1 74 

00 910 1=1, NS 75 

LL= LH I GH ( I ) 76 

910 XM< I )= X{LL ) 77 

DEBUG NB.NS, (XM( I ), 1=1, NS) 78 

DEBUG (LLOW(I), 1=1, NS) 79 

DEBUG (LHIGHt I) , 1=1, NS) 80 

C 81 

RETURN 82 

END 83 


IBFTC FFHGH 

1 

PROGRAM VARIABLES 2 

***************** 3 

4 

NSTRL - TRIAL NUMBER OF SEGMENTS 5 

NS - NEW SEGMENT COUNTER 6 

NPLFT - NUMBER OF POINTS LEFT THAT HAVE NOT BEEN ASSIGNEO 7 

TO A SUBSET 8 

N - INDEX OF HIGHEST POINT USED FOR LAGRANGE POLYNOMIAL 9 

N1 - INDEX OF SECOND HIGHEST POINT 10 

N2 - INOEX OF THIRD HIGHEST POINT - LOWEST POINT FOR QUADRATIC 11 

N 3 - INDEX OF LOWEST POINT FOR CUBIC 12 

NM - INDEX OF FIRST POINT TO BE TESTED 13 

******* 14 

A . 15 

R *- INTERMEDIATE VALUES TO SIMPLIFY CODING OF THE 16 

C * LAGRANGE POLYNOMIAL 17 

0 * 18 

******* 19 

YJ - Y AT X(J) EVALUATED BY LAGRANGE POLYNOMIAL 20 

21 

SUBROUTINE DHIX , Y,XM,LLOW,LHIGH, IX, IXM, IL) 22 

DIMENSION X t I X) , Y ( IX ) , XM ( I XM ) , LLOW I I L ) , LH I GH I I L ) 23 

RETURN 24 

25 

DETERMINE SUBSETS BY FORCE FITTING STARTING AT HIGH END OF DATA 26 

27 

ENTRY FFHIGHIX, Y , XM , L LO W , LH I GH , NX , NS , NB , M , TOL ) 28 

C 29 

NSTRL = NS 30 

NN = NS 31 

NS = 0 32 

800 NPLFT = NX 33 

DEBUG NPLFT, NSTRL 34 

DO 810 1=1, NSTRL 35 

LHIGHI I )= 0 36 

810 LLOW ( I )= 0 37 

LHIGH(NN)= NX 38 

DO 860 1=1, NN 39 

IF (NPLFT. LE.M) GO TO 870 40 

820 II = NSTRL-I+1 41 

NS = NS+1 42 

DEBUG I, II, NS, NPLFT 43 

N = NPLFT 44 

N 1 = N— 1 45 

N 2 = N-2 46 

N3 = N-3 47 

NM = N3-M+2 48 

NPLFT = NPLFT-M 49 

DEBUG I ,N,N1,N2,N3,NM 50 

DO 850 J= 1 , NM 51 

JJ = NM-J+1 52 

IF (J.EQ.l) LLOW( II ) = JJ 53 

DEBUG J, JJ.Xf JJ ) . Y< JJ ) 54 

IF (M.EQ.3) GO TO 830 55 
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A = (X(Nl)-X(JJ) >/( X( N)— X(N2) > 56 

B = (X(N)-X(JJ) )/(X(Nl)-X(N2) ) 57 

C = t X ( N2 ) — X ( JJ ))/ ( X< N)-X(Nl) ) 58 

YJ = Y(N2)*A*B-Y(N1 )*B*C+Y(N)»C*A 59 

DEBUG A,B,C,YJ 60 

GO TO 840 61 

830 A = (XIN2I-XI JJ) >/(X(N)-X(N3> > 62 

B = (XINl >-X< JJ I )/ I XI N21-XIN3) ) 63 

C = t X ( N ) —X (JJ) )/(X(Nl) —X I N3 ) ) 64 

D = ( X I N3 ) — X ( J J) ) / ( X (N 1 ) -X ( N2 ) ) 65 

YJ = Y(N3)*A*B*C-Y(N2)*D*B*(X(N)-X( JJ )) / ( X I N ) -X I N2 > >+YINl)*D*C. 66 

1 (X(N2)-X( JJ) )/ (X(N)-XI Nl) )-Y(N)»A* IXIN3) -X(JJ))»(X(N1)-X(JJ))/ 67 

2 (X(N)-XIN2> )/(X(N)-X(Nl) ) 68 

DEBUG AfBfCfDtYJ 69 

840 IF (YIJJ).EQ.O. ) GO TO 841 70 

IF ( A8S( l.-YJ/Y ( JJ) J.GT.TOL) GO TO 859 71 

LLOWUI ) = JJ 72 

GO TO 850 73 

841 IF ( A8S ( YJ-YI JJ )) .GT.TOL ) GO TO 859 74 

LLOWUI )= JJ 75 

850 NPLFT = NPLFT-1 76 

859 IF III.NE.l) LHIGFI I I — 1 )=LLQW( I I ) 77 

860 CONTINUE 78 

C 79 

870 L LOW ( II ) = 1 80 

DEBUG ( LLQWt I ), 1=1, NSTRL ) 81 

DEBUG (LHIGHI I) , 1=1, NSTRL) 82 

DO 880 I = I I , NSTRL 83 

IJ = 1-11*1 84 

LLOW ( IJ )= L LOW I I I 85 

LHIGHI I J ) = LHIGH (II 86 

JJ= LHIGH(IJ) 87 

XM ( I J ) = X ( J J ) 88 

880 CONTINUE 89 

N 8= NS — 1 90 

DEBUG ( XM( I I , 1=1, NS I 91 

DFBUGILLOWI I I ,1 = 1,NS) 92 

DEBUG (LHIGHI!) , 1 = 1 , N S ) 93 

C 94 

RETURN 95 

END 96 


SIBFTC REFT 

C 1 

C PROGRAM VARIABLES 2 

C ................. 3 

C 4 

C NST - NUMBER OF THE SUBSET FROM WHICH POINTS ARE BEING 5 

C CHECKED 6 

C NSS - NEW SUBSET COUNTER 7 

C 8 

C AA . 9 

C BB *- COEFFICIENTS FOR TESTING POLYNOMIAL, 10 

C CC • Y = AA + BB*X + CC*X**2 + DD*X*»3 11 

C OD * ( IF M = 2 , CD = 0 • I 12 

C ....... 13 

C I - INDEX OF POINT BEING TESTED 14 

C YI - Y EVALUATED AT XII) BY TESTING POLYNOMIAL 15 

C YJ - Y EVALUATED BY COE FF I E NTS FOR SEGMENT NST 16 

C NS - NUMBER OF SEGMENTS IN FIRST FIT 17 

C 18 

SUBROUTINE RF T( X. A , XM ,L LOW , LH l GH , I X , I A , I XM , I Ml ) 19 

DIMENSION X(IX), A ( I A . I Ml I , XM (I XM) .LLOW II A) , LH I GH II A I 20 

DOUBLE PRECISION A , AA , BB.CC.DD 21 

GO TO 140 22 

23 

CHECK IF DATA SHOULD BE REFITTED AND DETERMINE NEW SUBSETS 24 

25 

ENTRY REFIT (X,A,XM,LLOW,LH I'GH. NX.NS .NSS.M.TOL) 26 

DEBUG I LLOW (I), I-l.NS) 27 
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DEBUG (LHIGH(I). 1=1, NS) 28 

NST=2 29 

NSS= 1 30 

I S T= LH IGH ( 1 ) +1 31 

A A= ACl.l) 32 

BB= A ( 1 , 2 ) 33 

CC= A(l, 3) 34 

DD= 0 . DO 35 

I F ( M . EQ • 3 ) DD=A( 1, A) ' 36 

DEBUG AA,BB,CC,DD,NSS,NST 37 

C 38 

100 DO 130 1= 1ST, NX 39 

I F t X ( I) .GT.XM(NSTJ) NST= NST+1 AO 

Y I = AA+XC I)*(BB+X( I )*(CC+DD*X( I ))> A1 

Y J= A ( NST , M+ 1 ) A2 

DO 110 J=1»M A3 

IJ = M-J + l AA 

YJ= YJ»X( I )+A(NST, I JI A 5 

110 CONTINUE A6 

DEBUG NSS, NST,l.YI,YJ A7 

IF (YI.NE.O. ) GO TO 120 A8 

IF (ARS(YI-YJ).LE.A8S(TOL)) GO TO 130 A9 

GO TO 125 50 

120 IF (ABSI1.-YJ/YI I.LE.ABSITOLI) GO TO 130 51 

52 

END OF NSS SUBSET 53 

5A 

125 LHIGH(NSS)= 1-1 55 

IF (NSS.NE.NS) LL0W(NSS+1)=I-1 56 

I ST= I +M 57 

IF (IST.GT.NX) GO TO 135 58 

NSS= NSS+1 59 

A A= AINST.l) 60 

BB= AINST.2) 61 

CC = A ( N ST , 3 ) 62 

IFIM.FQ.3) DD= A ( NS T, A ) 63 

DEBUG I , NSS. NST , AA, BB,CC*DD 6A 

GO TO 100 65 

130 TOL = -ABS(TOL) 66 

135 LHIGH ( NSS ) = NX 67 

DEBUG (LLOWI I ), 1=1, NSS) 68 

DEBUG (LHIGH(I) , 1=1, NSS) 69 

70 

1 AO RETURN 71 

END 72 


IBFTC MINV 

1 

PROGRAM VARIABLES 2 

3 

A 

AIN - MATRIX TO EE INVERTED 5 

NN - ORDER OF AIN 6 

DET - VALUES OF THE DETERMINANT 7 

ERR - MAXIMUM DEVIATION OF ELEMENTS OF A I N« A I N ( I N VERSE ) 8 

9 

A - WORKING MATRIX 10 

N - NUMBER CF ROWS IN WORKING MATRIX 11 

JND - NUMBER OF COLUMNS IN WORKING MATRIX 12 

AK - VALUE OF THE FIRST ELEMENT IN PIVOTAL ROW 13 

FROM UNIT MATRIX 1A 

ERR 1 - SCALAR PRODUCT OF I ( TH I ROW OF AIN AND J(TH) COLUMN 15 

OF AIN( INVERSE) . ALSO, THE DEVIATION OF THE I.JITH) 16 

ELEMENT FROM UNIT MATRIX 17 

18 

MATRIX INVERSION BY GAUSSIAN ELIMINATION 19 

20 

SUBROUTINE M I NVRT ( A IN, NN, DE T, ERR ) 21 

DOUBLE PRECISION A I N, A , AK , DET , ERR ,E RR1 22 

DIMENSION A I N 1 A » A ) »-A! A , 81 23 
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24 

TRANSFER INPUT MATRIX (AIN) TO WORKING ARRAY (A) AND FILL 25 

REMAINDER OF WORKING ARRAY WITH UNIT MATRIX 26 

27 

N= NN 28 

JND= 2»N 29 

DO 110 1=1, N 30 

DO 100 J=1,N 31 

A ( I , J ) = AIN(I,J) 32 

JN= J + N 33 

A ( I , JN ) = 0 • DO 34 

IF(l.EQ.J) All , JN)= 1.00 35 

100 CONTINUE 36 

DEBUG ( A ( I , J ) , J=1 , JND ) 37 

110 CONTINUE 38 

39 

CALCULATE DETERMINANT AND ELIMINATE THE .BELOW THE DIAGONAL. 40 

ELEMENTS OF LEFT SIDE OF A 41 

42 

DET = 1 . DO 43 

DO 230 I=1,N 44 

DET = DET *A( I , I ) 45 

IF ( A ( I , I ) • EQ .0 • DO ) GO TO 600 46 

J ST = I 47 

AK = A( I, I ) 48 

DO 200 J= J ST, JND 49 

200 A ( I , J)= A( I , J )/ AK 50 

DEBUG I ,(A( I , J ) ,J=1,JND) 51 

IF(I.EQ.N) GO TO 3C0 52 

K ST= 1 + 1 53 

DO 220 K=KST ,N 54 

AK= A ( K , I ) 55 

DO 210 J= J ST, JND 56 

210 A ( K , J ) = A ( K , J ) — A ( I , J ) * AK 57 

DEBUG K, ( A1K,J ) , J= 1, JND) 58 

220 CONTINUE 59 

230 CONTINUE 60 

OEBUG DET 61 

62 

ELIMINATE THE ABOVE THE 0 1 AGONAL ELEMENTS OF LEFT SIDE OF A 63 

64 

300 DO 330 1=2, N 65 

K ND= 1-1 66 

DO 320 K= 1 , KND 67 

AK= A ( K , I I 68 

J ST= I 69 

DO 310 J= JST.JND 70 

A ( K, J ) = A ( K • J ) - A ( I , J ) • AK 71 

310 CONTINUE 72 

DEBUG K,(A(K,J) ,J=1,JND) 73 

320 CONTINUE 74 

330 CONTINUE 75 

DEBUG N, IAIN, J) ,J=1 , JND) 76 

77 

INVERSE IS IN RIGHT SIDE OF A. MULTIPLY INPUT MATRIX BY 78 

THE INVERSE AND FIND THE MAXIMUM DEVIATION OF ELEMENTS 79 

OF THE PROOUCT MATRIX FROM THE UNIT MATRIX 80 

81 

ERR= 0 . DO 82 

DO 420 1=1, N 83 

DO 410 J= 1 , N 84 

JN= J+N 85 

ERR1= 0.00 86 

DO 400 K=1,N 87 

400 ERR1= ERR1+AINI I ,K) «A(K, JN) 88 

ERR1= DABS ( ERR1 ) 89 

IFIJ.EQ.I) ERR1= ERR 1-1. DO 90 

IF(ERRl.GT.ERR) ERR= ERR1 91 

DEBUG I , J , ERR > ERR 1 92 

410 CONTINUE 93 

420 CONTINUE 94 

95 

TRANSFER INVERSE TO INPUT ARRAY 96 
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C 97 

00 5 10 1=1, N 98 

DO 500 J=1,N 99 

JN= J+N 100 

A INI I,J)= A( I , J N ) 101 

500 CONTINUE 102 

510 CONTINUE 103 

RETURN 104 

C 105 

C FOR SINGULAR MATRIX, RETURN NULL MATRIX 106 

C 107 

600 DO 620 1=1, N 108 

DO 610 J = 1 ,N 109 

610 A ( I » J ) = 0 . 00 110 

620 CONTINUE 111 

RETURN 112 

C 113 

END 114 


$IBFTC MATDEF 

C 1 

C PROGRAM VARIABLES 2 

c 3 

C 4 

C T - ONE BLOCK OF MATRIX PRODUCT X-TRANSPOSE »W*X 5 

C TT - POWER OF XII) IN FORMATION OF T 6 

C DET - DETERMINANT OF T 7 

C ERR - VALUE THAT MEASURES ACCURACY OF T-INVERSE 8 

C 9 

SUBROUTINE DEF I XX , Y Y, W , XM , LLOW , LHI GH, XWX, YWX , C , I X , I XM , I L , I M 1 , I M, 10 

1 INI) 11 

DIMENSION XXI IX ) ,YY (IX) ,W( IX) ,XM( IXM) ,LLUW( IL I.LHIGHI IL) , 12 

1 YWX( IL, IM1 ) ,XWXI IL, IM1, IM1) ,C( INI , IM, IM1) ,T(4,4) 13 

DOUBLE PRECISION YW X, XW X , C , T , TT , DET , ERR 14 

GO TO 400 15 

C 16 

C DEFINE THREE MULTIDIMENSIONAL MATRICES REQUIRED FOR SOLUTION 17 

C OF VECTOR A-TRANSPOSE 18 

C 19 

ENTRY DEFMATI XX, YY, W, XM , LLOW, LH I GH, NX , N S, M, XWX, YWX, C > 20 

C 21 

C DEFINE MATRIX OF CONSTRAINTS IF THERE IS MORE THAN ONE SEGMENT 22 

C 23 

100 NN = NS-1 24 

Ml = M+l 25 

MM = M 26 

IF (NS.EQ.l) GO TO 200 27 

DO 120 N= 1 , NN 28 

C 29 

C FIRST ROW 30 

C 31 

C ( N, 1 , 1 ) = l.DO 32 

DO 110 K=2, Ml 33 

110 C(N.l.K) = C1N, 1,K-1)*XM(N) 34 

C 35 

C SECCND ROW 36 

C 37 

C IN, 2, l ) = O.DO 38 

C ( N, 2 , 2 ) = l.DO 39 

CIN.2,3) = 2 . DO * XM ( N) 40 

IF (MM.EQ.2I GO TO 120 41 

CIN.2.4) = 3.D0*XM( N)»*2 42 

C 43 

C THIRD ROW 44 

C 45 

CIN.3,1) = O.DO 46 

C ( N, 3, 2 ) = 0.00 47 

C <N,3,3) = 2. DO 48 

C I N , 3 ,4 ) = 6.D0*XM(N), 49 
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120 CONTINUE 50 

DO 130 N=1,NN 51 

DE8UG N 52 

DO 130 J = 1 , MM 53 

DEBUG J, (CtN,J,K),K=l, Mil 54 

130 CONTINUE 55 

C 56 

C DEFINE MATRIX XWX AND VECTOR YWX 57 

C 58 

200 NN = NS 59 

DO 300 N=1,NN 60 

KST= LLOWIN) 61 

KND- LHIGH(N) 62 

DEBUG KST.KND 63 

DO 210 J= 1 , Ml 64 

T( 1, J) = 0.00 65 

T(J,M1) = O.DO 66 

210 YWXIN.J) = O.DO 67 

DO 240 K=KST,KND 68 

TT = W(K> 69 

DO 220 J- 1 1 Ml 70 

TI1,J) = T ( 1 , J ) +T T 71 

YWXIN.J) = YWX(N,J)+TT*YY(K) 72 

220 TT = TT*XX(K) 73 

DO 230 1=2, Ml 74 

T ( I, Ml) = T ( I ,M 1 ) +TT 75 

TT = TT*XX(K) 76 

230 CONTINUE 77 

240 CONTINUE 78 

C 79 

DO 270 1=2, Ml 80 

DO 260 J= 1 , MM 81 

260 T I I , J) = T ( 1-1, J + l) 82 

270 CONTINUE 83 

DEBUG N 84 

DO 271 1=1, Ml 85 

DEBUG I , ( T ( I , J ) , J= 1 »M 1 ) 86 

271 CONTINUE 87 

DEBUG ( YWXIN, I) , 1 = 1, Ml ) 88 

C 89 

CALL MINVRT(T,M1,0ET,ERR> 90 

DEBUG GET, ERR 91 

DO 290 1=1, Ml 92 

0E81JG I.ITI I, J) , J = 1 »M 1 ) 93 

DO 290 J= 1 , M 1 94 

XWXI N, I , J 1 = T( I , J ) 95 

290 CONTINUE 96 

300 CONTINUE 97 

C 98 

400 RETURN 99 

END 100 


tIBFTC ASL V 
C 

C PROGRAM VARIABLES 

C **».*****»**»»*». 

C B - MATRIX PRODUCT C» ( X-TR ANSPOSE • W»X ) *C-TRANS POS E 

C N - ROW INDEX OF SUBMATRICES OF B 

C L - COLUMN INDEX OF SUBM ATR ICES OF B (L=l,2,3) 

C I - ROW INDEX OF ELEMENTS OF THE SUBMATRICES BIN.L) 

C J - COLUMN INDEX OF ELEMENTS OF THE SUBMATRICES BIN.L) 

C SIGN - PLUS OR MINUS 1 - CHANGES SIGN OF MATRIX PRODUCT 

C TT - INTERMEDIATE MATRIX IN THE DOUBLE MULTIPLICATION 

C T - INTERMEDIATE MATRIX IN THE DOUBLE MULTIPLICATION 

C BINV - INVERSE OF B 

C E - MATRIX E FROM SOLUTION FOR B-INVERSE IN APPENDIX B 

C DELTA - MATRIX DELTA FROM SOLUTION FOR B-INVERSE 

C TS - INTERMEDIATE VALUE RELATED TO THE IDENTITY SUBMATRICES 

C OF APPENDIX B 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 
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V - VECTOR V OF APPENDIX B 17 

VV - VECTOR VV OF APPENDIX B 18 

19 

SUBROUTINE SL VI C , XWX, YWX, A , B , B I NV, IC , I M , I Ml , IL) 20 

DIMENSION CIIC, I M, I Ml) , XWX ( I L . I Ml , ! Ml ) . YWX I I L » I Ml ) , A 1 1 L , I Ml ) , 21 

1 B(IC,IM,IM,IM) f BINV(IC,IC,IM f IM) 22 

DIMENSION Tf4,4),TT<4,4),E(4,4),DELTA(3,3),V(4),VV(4) 23 

DOUBLE PRECISION C, XWX, YWX, A , B , B I NV, T , TT, E , DELTA , V, VV, TS, DE T ,D I V , 24 

1 SIGN 25 

GO TO 940 26 

27 

SOLVE MATRIX EQUATION FOR A-TRANSPOSE 28 

29 

ENTRY ASOLVECC,XWX,YWX,A,NS,M) 30 

31 

FOR ONE SEGMENT, DO A SIMPLE LEAST SQUARES FIT 32 

33 

IF (NS.GT.lt GO TO 100 34 

Ml = M+l 35 

DO 95 1=1, Ml 36 

All, I) = O.DO 37 

DO 90 K = 1 , M 1 38 

All, I) = A( 1,1 t+XWX (1, I ,K)*YWX( 1»KI 39 

90 CONTINUE 40 

95 CONTINUE 41 

RETURN 42 

43 

DEFINE B MATRIX 44 

45 

100 NN = NS-1 46 

MM = M 47 

Ml = M+l 48 

DO 200 N= 1 » NN 49 

DO 190 L= 1 , 3 50 

NIND = (L-1I/2+N 51 

GO TO ( 1,2.3) ,L 52 

1 IF (N.EQ.l) GO TO 190 53 

JIND = N-l 54 

GO TO 120 55 

C 56 

2 SIGN = 1 . DO 57 

JIND = N 58 

DO 110 1=1, Ml 59 

DO 110 J= 1 , Ml 60 

T ( I , J ) = XWXIN, I , J) +XWXIN+1, I , J) 61 

110 CONTINUE 62 

GO TO 140 63 

C 64 

3 IF (N.EO.NN) GO TO 200 65 

JIND = N+l 66 

120 SIGN = -l.DO 67 

DO 130 1=1, Ml 68 

DO 130 J= 1 , M 1 69 

T ( I , J ) = XWXININO, I ,J ) 70 

130 CONTINUE 71 

C 72 

140 DO 160 1=1, Ml 73 

DO 150 J= 1 , MM 74 

TT II , J) = O.DO 75 

DO 150 K= 1 , M 1 76 

150 TT ( I , J ) = TT( I, J)+T(I,KI*C(JIND,J,K) 77 

DEBUG I , I TT 1 1 , J ) , J= 1, MM ) 78 

160 CONTINUE 79 

C 80 

DO 180 1=1, MM 81 

DO 175 J= 1 , MM 82 

B(N,L,I,J) = O.DO 83 

DO 170 K=1 , Ml 84 

170 B(N,L,I,J) = B(N,L, I, J I +C I N, I ,K ) *TT ( K , J ) 85 

175 BCN.L.I.J) = BIN.L, I, J)*SIGN 86 

180 CONTINUE 87 

190 CONTINUE 88 

200 CONTINUE 89 

C 90 
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DO 370 K=1,MM 164 

DELTA ( I f J ) = DELTA! I ,J )+8(LL,l, 1 ,K)*TIK, J) 165 

370 E ( I • J ) = £{ I,J)+B(LL,1,I,K)*TT(K.J) 166 

E { 1 » J ) = 6ILL.2, I , JI-EI I , J) 167 

TS= 0.D0 168 

IFILL.EQ.N.AND. I.EQ.J I TS=1.D0 169 

DELTA(I,J)= TS— DELT A( I » J ) 170 

380 CONTINUE 171 

390 CALL MINVRTI E ,MM f DE T, DE V ) 172 

396 CONTINUE 173 

C 174 

400 DO 420 1=1, MM 175 

DO 410 J= 1 , MM 176 

DO 410 K= 1 , MM 177 

DELT A I l » J > = DELTA ( I ,J >-B(L,3, I,K) *BI NV t L + l , N , K , J > 178 

410 CONTINUE 179 

420 CONTINUE 180 

C 181 

DO 440 1=1, MM 182 

DO 430 J=1 ,MM 183 

BINVIL.N, I, J)= 0.00 184 

DO 430 K=1 ,MM 185 

BINVIL.N, I, J)= BINVIL.N, I, JI+E(I,K) •DELTAIK.J) 186 

430 CONTINUE 187 

DEBUG L.N. (BINVIL.N, I, J),J=1, MM) 188 

440 CONTINUE 189 

C 190 

IF (L.EQ.l) GO TO 450 191 

L = L-l 192 

GO TO 310 193 

450 CONTINUE 194 

195 

CARRY GUT MATRIX MULTIPLICATIONS FOR A-TRANSPOSE 196 

197 

500 DO 930 N= 1 , NS 198 

DO 510 L= 1 , M 1 199 

VIL) = YWXIN.LI 200 

510 CONTINUE 201 

DEBUG N,(V(L),L=1,M1> 202 

DO 850 J J = 1 , NS 203 

204 

FORM MATRIX OIJJ.N) BY 4 SEPARATE TECHNIUUES 205 

206 

IFIN.NE.l) GO TO 550 207 

NDUM= 1 208 

IFIJJ.NE.l) GO TO 530 209 

JDUM= 1 210 

SIGN = 1 • DO 211 

GO TO 570 212 

530 IFIJJ.NE.NS) GO TO 540 213 

JDUM= NS— 1 214 

SIGN = — 1 • DO 215 

GO TO 570 216 

540 S I GN= 1 . DO 217 

GO TO 630 218 

550 IFIN.NE.NS) GO TO 561 219 

NDUM= NS— 1 220 

IFIJJ.NE.l) GO TO 560 221 

J DUM= 1 222 

SIGN= - 1 . DO 223 

GO TO 570 224 

560 IFIJJ.NE.NS) GO TO 563 225 

J DUM= NS— 1 226 

SIGN = 1 . DO 227 

GO TO 570 228 

561 IF IJJ.NE.l) GO TO 562 229 

JDUM = 1 230 

SIGN = 1 . DO 231 

GO TO 615 232 

562 IF (JJ.NE.NS) GO TO 680 233 

JDUM = NS— 1 234 

SIGN = — 1 . DO 235 

GO TO 615 236 
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FIND INVERSE OF B 91 

92 

DO 450 N= 1 y NN 93 

DO 280 L=1,NN 94 

IF IL.NE.l) GO TO 220 95 

DO 210 1=1, MM 96 

DO 210 J= 1 , MM 97 

E ( I , J ) = B ( 1, 2, I,J) 98 

DELT A ( I , J ) = 0. DO 99 

IF (L.EQ.N.AND. I.EQ.J) DELTA ( I , J ) =1 . DO 100 

210 CONTINUE 101 

GO TO 270 102 

C 103 

220 DO 240 1=1, MM 104 

DO 230 J= 1 , MM 105 

T ( I, J ) = 0 . DO 106 

TT(I.J) = 0 • DO 107 

DO 230 K= 1 , MM 108 

TII.J) = T< I,J)+E(I,K)*DELTA(K,JJ 109 

TT(I,J) = TT( I, J )+E(I ,K )*B(L-1,3,K, J) 110 

230 CONTINUE 111 

240 CONTINUE 112 

C 113 

DO 265 1=1, MM 114 

DO 260 J=1 , MM 115 

E(I,J)=0.D0 116 

DELT A ( I , J ) = 0. DO 117 

DO 250 K = 1 * MM 118 

E(I,J> = E(I,J)+B(L,1,I,K)*TT(K,J) 119 

250 DELTA ( I , J ) = DELTA! I,J)+B(L,1,I,K)*T(K,J) 120 

TS = 0 . DO 121 

IF (L.EQ.N.AND. I.EQ.J) TS=1.D0 122 

DELT A ( I , J ) = TS— DELTA ( I • J ) 123 

E ( I , J ) = B ( L , 2 , I , J ) — E ( I , J ) 124 

260 CONTINUE 125 

265 CONTINUE 126 

C 127 

270 CALL MINVRT(E,MM,OET,DEV) 128 

280 CONTINUE 129 

C 130 

DO 300 1=1, MM 131 

DO 290 J= 1 , MM 132 

BINVINN.N, I,J) = 0.00 133 

DO 290 K=1,MM 134 

290 BINVINN.N, I, J)= B IN V( NN , N , I , J ) + E ( I ,K ) *DELTA ( K , J ) 135 

DEBUG NN,N, (BINVINN.N, I »J) ,J=1.MM( 136 

300 CONTINUE 137 

C 138 

IF (NN.'EQ.l) GO TO 500 139 

L = NN-1 140 

310 LAST = L 141 

DO 396 LL= 1 , L AS T 142 

320 IF(LL.NE.l) GO TO 340 143 

DO 330 1=1, MM 144 

DO 330 J= 1 , MM 145 

E ( I , J ) = B ( 1 , 2, I , J ) 146 

DELTAU ,J>= O.DO 147 

IFILL.EQ.N. AND. I.EQ.J) DELTA(I,J)= l.DO 148 

330 CONTINUE 149 

GO TO 390 150 

C 151 

340 DO 360 1=1, MM 152 

DO 350 J= 1 , MM 153 

T ( I , J ) = O.DO 154 

TT ( I , J ) = O.DO 155 

DO 350 K= 1 , MM 156 

T ( I , J ) = T( I, J )+E ( I,K)*DELTA(K, J) 157 

350 T T ( I , J ) = TTU, J)+E(I,K)*B(LL-1,3,K,J) 158 

360 CONTINUE 159 

DO 380 1=1, MM 160 

DO 380 J=1 , MM 161 

DELT A ( I , J)= O.DO 162 

E( I, J)= O.DO 163 
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563 SIGN = -1.00 237 

GO TO 630 238 

C 239 

570 DO 590 1=1, Ml 240 

DO 580 J=1,MM 241 

T ( I , J I = 0.00 242 

DO 580 K=1,MM 243 

T I I , J ) = T( l , J )+C I JDUM.K, I )*BINV! JDUM.NDUM ,K, J) 244 

580 CONTINUE 245 

DEBUG ITU.JI, J = 1 , KM ) 246 

590 CONTINUE 247 

C 248 

DO 610 1=1, Ml 249 

DO 605 J=1 ♦ Ml 250 

TT ( I , J)= 0 . DO 251 

DO 600 K=1.MM 252 

TT ( I , J ) = TT ( I , J ) +T(I ,K)*C(NDUM,K,J) 253 

600 CONTINUE 254 

TT ( I , J ) = TT(I,J)«SIGN 255 

605 CONTINUE 256 

DEBUG ( TT ( I » J) » J= 1, Ml ) 257 

610 CONTINUE 258 

GO TO 800 259 

C 260 

615 DO 617 1=1, MM 261 

DO 616 J= 1 , Ml 262 

T ( I , J ) = 0 • DO 263 

DO 616 K= 1 , MM 264 

T l I , J ) = Tt I , J) +BINVI JOUM.N, I,K l*C(N,K, JJ-BINVI JDUM.N-l, I ,K )» 265 

1 C (N— 1 , K, J ) 266 

616 CONTINUE 267 

DEBUG (T( I , J) ,J=1,M1) 268 

617 CONTINUE 269 

C 270 

DO 620 1=1, Ml 271 

DO 619 J = 1 , Ml 272 

TTI I , J) = 0.00 273 

DO 618 K= 1 , MM 274 

TTII.J) = TTI I, J )+C (JDUM.K, I )*T1K, J I 275 

618 CONTINUE 276 

TTII.J) = TTI I, J )*SIGN 277 

619 CONTINUE 278 

DEBUG (TTI I , J), J=l, Ml ) 279 

620 CONTINUE 280 

GO TO 800 281 

C 282 

630 DO 650 1=1, Ml 283 

DO 640 J = 1 , MM 284 

T I I , J ) = O.DO 285 

DO 640 K= 1 , MM 286 

TII.J) = T I I , J) +C I J J.K, I I *BINV( JJ.NDUM.K, JI-C I JJ-1, K, I >• 287 

1 BINV I JJ-1 , NDUM.K, J ) 288 

640 CONTINUE 289 

DEBUG (TII.J), J = 1 , MM ) 290 

650 CONTINUE 291 

C 292 

DO 670 1=1, Ml 293 

DO 665 J = 1 , M 1 294 

TT I I » J ) = O.DO 295 

DO 660 K= 1 , MM 296 

TT I I , J) = TTI I,J )+T< I,K) *C(NDUM,K, J) 297 

660 CONTINUE 298 

TTII.J) = TTI I, J )»S IGN 299 

665 CONTINUE 300 

DEBUG (TTII.J), J = 1 ,M 1 ) 301 

670 CONTINUE 302 

GO TO 800 303 

C 304 

680 DO 700 1=1, Ml 305 

DO 690 J=1 , MM 306 

T I I , J ) = O.DO 307 

DO 690 K=1 , MM 308 

TII.J) = T( I, Jl+CI JJ-l.K. I )*8INVt JJ-l.N-l.K, Jl-CI JJ.K, I )» 309 

1 BINVI JJ,N-1,K,J) 310 
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690 CONTINUE 311 

DEBUG <T( 1 1 J ) ,J = 1 ,MM) 312 

700 CONTINUE 313 

C 314 

DO 720 1=1, Ml 315 

DO 710 J= 1 , Ml 316 

TT ( I , J ) = 0.00 317 

DO 710 K=1,MM 318 

TT 1 1 , J ) = TTII.JJ + TI I,K)*C(N-1,K,J) 319 

710 CONTINUE 320 

DEBUG ITT1I.J). J=1.M1) 321 

720 CONTINUE 322 

C 323 

DO 740 1=1, Ml 324 

DO 730 J = 1 , MM 325 

T< I, J)= 0.00 326 

DO 730 K= 1 , MM 327 

TII.J) = T< I , J) +CI JJ, K, I )»BINV( JJ.N.K, J )-C( JJ-l.K.I >* 328 

1 BINVI JJ-l.N.K.J) 329 

730 CONTINUE 330 

DEBUG <T( I,J),J=1,MM) 331 

740 CONTINUE 332 

C 333 

DO 760 1 = 1, Ml 334 

DO 750 J= 1 , Ml 335 

DO 750 K= 1 , MM 336 

TT ( I , J ) = TT l I , J ) + T( I ,K)*C(N,K,J) 337 

750 CONTINUE 338 

DEBUG ( TTI I ,JI, J=1,M1 ) 339 

760 CONTINUE 340 

341 

MATRIX D(JJ.N) IS STORED IN TT 342 

343 

800 DO 820 L= 1 , Ml 344 

V V ( L ) = 0 . DO 345 

DO 810 K= 1 , Ml 346 

VV(L) = VV(L)+YWX(JJ,K)*XWX(JJ,K,L) 347 

810 CONTINUE 348 

820 CONTINUE 349 

DEBUG JJ, IVV(L) *L = 1 ,M1) 350 

DO 840 L=1 , Ml 351 

DO 830 K= 1 , M 1 352 

V t L ) = V(L)-VVIK)*TT(K,U 353 

830 CONTINUE 354 

840 CONTINUE 355 

DEBUG JJ, (V(U,L = 1, Ml ) 356 

850 CONTINUE 357 

358 

FINAL MULTIPLICATION FOR A-TRANSPOSE 359 

360 

900 DO 920 1=1, Ml 361 

AIN, I ) = 0.00 362 

DO 910 K=1 , Ml 363 

AIN. I) = A(N,I)»V(K)*XWX(N,K, I) 364 

910 CONTINUE 365 

920 CONTINUE 366 

DEBUG N 367 

DEBUG (AIN, I), 1=1, Ml) 368 

930 CONTINUE 369 

C 370 

940 RETURN 371 

END 372 



Subroutine TRANSF 



Subroutine BTRANS 






Subroutine ORDER 
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Subroutine DIVXM 
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Subroutine DIVNS 
































Subroutine REFIT 






Subroutine MINVRT 



Transfer input array to left side of working array and 
fill right side of working array with unit matrix 


Diagonalize left side of working array and calculate determinant) 


Eliminate above diagonal elements on left side of 
(inverse is now in right side of working 




ERR * maximum |(AA _1 )jj - Iyj 




Transfer inverse to input array 


working array 
array) 


Return] 









Subroutine ASOLVE 






















APPENDIX F 


COMPUTER INPUT AND OUTPUT SHEETS FOR SAMPLE PROBLEM 1 
IN PUT FOR SAMPLE PROBLEM 1 






Card column 





1-6 

7-12 

13-18 

19-24 

25-30 

31-36 

37-42 

43-48 

49-60 

61-66 


SAMPLE 

PROBL 

EM 1 







2 45 

0 2 

0 F 

F F 

F .r 

01 





(6F12 

.5) 









8 

. 86 

5 

.556 



10 

.98 

7 

.6798 

12 

.90 

13 

.116 



16 

.81 

11 

.96 7 

60 

. 805 

43 

.043 



120 

.86 

86 

.038 

180 

.88 

116 

.49 



240 

. 82 

161 

.03 

248 

.09 

158 

.00 



255 

. 89 

159 

.89 

265 

. 94 

168 

.00 



275 

.78 

174 

.92 

279 

.30 

181 

.84 



300 

.85 

210 

.51 

3 60 

.89 

244 

.47 



374 

. 74 

238 

.03 

406 

. 99 

260 

.51 



415 

.79 

275 

.91 

445 

.53 

287 

.06 



480 

.48 

320 

.67 

506 

.57 

324 

.27 



513 

.07 

327 

.02 

573 

.47 

372 

.89 



600 

. 79 

400 

*49 

626 

.50 

407 

.06 



644 

.26 

414 

.99 

699 

.31 

455 

.99 



722 

.16 

489 

.77 

900 

. 80 

598 

.51 



968 

.34 

624 

.53 

2865 

.4 

1837 

. 9 



3598 

.7 

2386 

.4 

3957 

. 1 

2554 

.0 



4281 

.0 

2821 

.0 

7200 

. 0 

4349 

.0 



8503 

.6 

5426 

.0 

11750 

.0 

7009 

.0 



12730 

.0 

7587 

.0 

15640 

. 0 

8788 

.0 



16940 

.0 

9427 

.5 

19550 

.0 

10360 

.0 



21160 

.0 

11060 

.0 

23440 

.0 

11720 

.0 



31230 

.0 

14020 

.0 

47300 

0 

14970 

.0 







( 12F6 

0 ) 









200 . 

0 

7000. 











OUTPUT FOR SAMPLE PROBLEM 1 
Linear Fit 


IND. VAR. 

oep.var. 

CALC. FUNC. 

DEVIATION 

RELATIVE ERROR 

8.85999990 

5.55599999 

473.876369 

-468.320366 

-0.98827541 

10.9799999 

7.67979997 

474.759743 

-467.079941 

-0.98382381 

12. 8999997 

13.1160001 

475.559776 

-462.443775 

-0.97241987 

16. 8099997 

11.9670000 

477.189011 

-465.222008 

-0.97492187 

60. 8049994 

43.0430002 

495.521046 

-452.478043 

-0.91313587 

120.859998 

86.038000 1 

520.545029 

-434.507027 

-0.83471554 

180.879997 

116.490COG 

545.554428 

-429.064426 

-0.78647409 

240.819998 

161.030001 

570.530502 

-409.500500 

-0.71775391 

248.089996 

158.000000 

573.559792 

-415.559792 

-0.72452741 

255.889997 

159.889999 

576.809937 

-416.919937 

-0.72280297 

265.939995 

168.000000 

580.997612 

-412.997612 

-0.71084218 

275.779995 

174.920000 

585.097786 

-410. 1 77784 

-0. 70104142 

279.299995 

181.840000 

586. 564514 

-404.724514 

-0.68999147 

300.849995 

210.510000 

595.544067 

-385.034065 

-0.64652489 

360.889996 

244.469999 

620.561806 

-376.091805 

-0.60605051 

374.739998 

238.03000 1 

626.332886 

-388.302883 

-0.61996247 

406.989994 

260.509998 

639. 770958 

-379.260960 

-0.59280740 

415.789997 

275.910000 

643.437782 

-367.527782 

-0.57119397 

445.529991 

287 .06000 1 

655.829979 

-368.769978 

-0. 56229509 

480.479996 

320.669998 

670.393097 

-349.723099 

-0.52166870 

506.569996 

324.270000 

681.264397 

-356.994396 

-0.52401739 

513.069992 

327.020000 

683.972847 

-356.952847 

-0.52188160 

573.469994 

372.889999 

709. 140587 

-336.250587 

-0.47416633 

600.789993 

400.490002 

720.524406 

-320.034405 

-0.44416872 

626.499992 

407.060001 

731.237366 

-324.177364 

-0.44332713 

644.259995 

414.990002 

738.637680 

-323.647678 

-0.43816838 

699.309990 

455.99G002 

761.576164 

-305.586163 

-0.40125489 

722.159996 

489.770000 

771.097404 

-281.327404 

-0.36484029 

900.799995 

598.510002 

845.533920 

-247.023918 

-0.29215140 

968.339981 

624.529999 

873.676781 

-249.146782 

-0.28517043 

2865.39996 

1837. 85999 

1664.15222 

173.747772 

0.10440617 

3598.69998 

2386.39999 

1969. 70695 

416.693039 

0.21155077 

3957.09998 

2554.00000 

2119.04666 

434.953339 

0.20525897 

4280.99994 

2821.0000C 

2254.01074 

566.989258 

0.25154683 

7199.99994 

4349. OCOOO 

3475.31274 

878.687256 

0.25320117 

8503.59985 

5426. OOCOO 

4013.50259 

1412.49741 

0.35193633 

11749.9999 

7009.00000 

5366.22705 

1642. 77295 

0.30613183 

12729.9999 

7587. OOCOC 

5774.57782 

1812.42218 

0.31386228 

15639.9998 

8788. OOOOC 

6987.12964 

1800. 87036 

0.25774108 

16939.9998 

9427. 50000 

7528.81946 

1898. 68054 

0.25218835 

19549.9998 

10360 .0000 

8616.36584 

1743. 63416 

0.20236306 

21159.9998 

11060.0000 

9287.22791 

1772. 77209 

0.19088280 

23439.9998 

11720.0000 

10237.2684 

1482. 73157 

0.14483663 

31229.9998 

14020. OOOC 

13483.2406 

536.759399 

0.39809376E-01 

47300.0000 

14970.0000 

20179.3601 

-5209.36011 

-0.25815289 


ThB regression equation for the above is 

Y = AO ♦ A1 X + ... 

THE PARAMETERS (AO-All ARE 

470.184547 0.41668448 

THE VARIANCE= 1346098.9 STANDARD DEVIATION= 1160.2150 

OETERMINANT= 86.72037 
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Parabolic Fit 


IND. VAR. 

DEP.VAR. 

CALC. FUNC. 

DEVIATION 

RELATIVE ERROR 

8.85999990 

5.55599999 

-7.82075906 

13.3767591 

-1.71041697 

10.9799999 

7.67979997 

-6.36010355 

14.0399035 

-2.20749605 

12. 8999997 

13.1160001 

-5.03730673 

18.1533067 

-3.60377234 

16.8099997 

11.9670000 

-2.34366494 

14. 3106649 

-6.10610527 

60.8049994 

43.0430002 

27.9484172 

15.0945830 

0.54008722 

120.859998 

86.038000 1 

69.2492838 

16.7887163 

0.24243884 

180.879997 

116.490C00 

110.469475 

6.02052498 

0. 54499444 E-0 1 

240.819998 

161.03000 I 

151.578245 

9.45175552 

0 • 62 355620 E-0 1 

248.089996 

158.000000 

156.560404 

1.43959618 

0.91951485E— 02 

255.889997 

159.889999 

161.904854 

-2.01485443 

-0. 12444682E-01 

265.939995 

168.000000 

168.789560 

-0.78956032 

-0.46777793E-02 

275.779995 

174.920000 

175.528872 

-0.60887146 

-0.34687824E-02 

279.299995 

181.840000 

177.939314 

3.90068626 

0.21921441E-01 

300.849995 

210.510C0C 

192.692173 

17.8178272 

0. 9246783 IE-0 1 

360.889996 

244.469999 

233.756330 

10.7136688 

0.45832636E-01 

374.739998 

238.030001 

243.220953 

-5. 19095230 

-0.21342537E-01 

406.989994 

260.509998 

265.247829 

-4.73783112 

-0. 17861903E-01 

415.789997 

275.910000 

271.255432 

4.65456772 

0.17159353E-01 

445.529991 

287 .06000 1 

291.549381 

-4.48937988 

-0. 15398352E-01 

480.479996 

320.669998 

315.380772 

5.28922653 

0.16770923E— 01 

506.569996 

324.2 7000C 

333.158276 

-8.88827515 

-0.26678836E-01 

513.069992 

327.02C0OC 

337.585655 

-10.5656548 

-0.31297700E-01 

573.469994 

372.889999 

378.694500 

-5.80450058 

-0. 15327660E-01 

600.789993 

400.490002 

397.269943 

3.22005844 

0.81054670E-02 

626.499992 

407.060001 

414.740002 

-7.68000031 

-0. 18517626E-01 

644.259995 

414.990002 

426.801937 

-11.8119354 

-0.27675449E-01 

699.309990 

455.990002 

464.158379 

-8.16837692 

-0.17598254E-01 

722.159996 

489. 7 7000C 

479.650211 

10.1197891 

0.21098269E-01 

900.799995 

598.510002 

600.481705 

-1.97170258 

-0.32 8 35348 E-0 2 

968.339981 

624.529999 

646.034935 

-21.5049362 

-0.33287575E-01 

2865.39996 

1837.89999 

1896.25816 

-58.3581696 

-0.30775435E-01 

3598.69998 

2386. 39999 

2364.37613 

22.0238647 

0 .93148735E— 02 

3957.09998 

2554.00000 

2590.09534 

-36.0953369 

-0. 139359UE-01 

4280.99994 

2821.0000C 

2792.35071 

28.6492920 

0. 10259919E-01 

7199.99994 

4349. OOOOC 

4540.73303 

-191.733032 

-0.42225128E-01 

8503.59985 

5426. OOCOO 

5278.30939 

147.690613 

0. 27 98 0666 E-0 1 

11749.9999 

7009.00000 

6999.10071 

9.89929199 

0. 14143663E-02 

12729.9999 

7587. 00000 

7486.02826 

100.971741 

0. 13488026E-01 

15639.9998 

8788.00000 

8842.99243 

-54.9924316 

-0.62187582E-02 

16939.9998 

9427. 50000 

9406.20850 

21.2915039 

0.22635586E-02 

19549.9998 

10360.0000 

10456.8165 

-96.8165283 

-0.92 587001 E-02 

21159.9998 

11060. COOO 

11051.5273 

8.47265625 

0.76665025E-03 

23439.9998 

11720.0000 

11824.0630 

-104.062988 

-0. 88009500E-02 

31229. 9998 

14020.0000 

13847.4071 

172.5 92*896 

0. 12463914E-Q1 

47300.0000 

14970. OCOO 

15009.6951 

-39.6950684 

-0.2 6 44 6 28 5 E-0 2 

THE REGRESSION 

EQUATION FOR THE ABOVE 

IS 



Y = AO + A1 X ♦ 

• • • 




THE PARAMETERS 

(A0-A2) ARE 




-13.9259609 

0.68914430 

— 0.78545492E— 05 




THE VARIANCE* 3248.4323 STANDARD DEVIATION* 56.995019 

DETERMINANT* 15.83930 
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Cubic Fit 


INO. VAR. 

8.85999990 

10.9799999 

12. 8999997 

16.8099997 

60.8049994 

120.859998 

180.879997 

240.819998 

248.089996 

255.889997 

265.939995 

275.779995 

279.299995 

300.849995 

360.889996 

374.739998 

406.989994 

415.789997 

445.529991 

480.479996 

506.569996 

513.069992 

573.469994 

600.789993 

626.499992 

644.259995 

699.309990 

722.159996 

900.799995 

968.339981 

2865. 39996 

3598.69998 

3957.09998 

4280.99994 

7199.99994 
8503. 59985 

11749.9999 

12729.9999 

15639.9998 

16939.9998 

19549.9998 

21159.9998 

23439.9998 

31229.9998 
47300.0000 


THB REGRESSION 

Y=A0+AIX+ 

THE PARAMETERS 
-7.98756391 

THE VARIANCE 
DETERMINANT* 0 


DEP.VAR. 

CALC. FUNC. 

DEVIATION 

RELATIVE ERROR 

5.55599999 

-1.94968595 

7.50568593 

-3.84968969 

7.67979997 

-0.50512648 

8.18492639 

-16.2037168 

13.1160001 

0.80309701 

12.3129030 

15.3317754 

11.9670000 

3.46707362 

8.49992633 

2.45161402 

43.0430002 

33.4265404 

9.61645985 

0.28768935 

86.038C001 

74.2768517 

11.7611485 

0.15834204 

116.490000 

115.050751 

1.43924904 

0. 12509688E-01 

161.030001 

155.717804 

5.31219673 

0.34114254E-01 

158.000000 

160.646660 

-2.64665985 

-0. 16475038E-01 

159.889999 

165.933983 

-6.04398346 

-0.36424024E-01 

168.000000 

172.745182 

-4.74518204 

-0.27469258E-01 

174.92000C 

179.412632 

-4.49263191 

-0.25040778E-01 

181. 84000G 

181.797392 

0.42608261E-01 

0.23437223E-03 

21O.51C0OC 

196.393326 

14.1166744 

0.71879603E— 01 

244.469999 

237.022972 

7.44702721 

0.31419010E-01 

238.0 3000 1 

246.387926 

-8.35792542 

-0.33921814E-01 

260.509998 

268.183537 

-7.67353821 

-0.28613010E-01 

275.910000 

274. 128231 

1. 781 76880 

0.64997639E-02 

287.06000 1 

294.210201 

-7.15019989 

-0.24303032E-01 

320.669998 

317. 793720 

2.87627792 

0. 90507701 E-02 

324.2 70C0C 

335.387054 

-11.1170540 

-0.33146938E-01 

327 .020000 

339.768669 

-12. 7486687 

-0.37521613E-01 

372.889999 

380.454430 

-7.56443024 

-0. 19882618E— 01 

400 .490002 

398.839809 

1.65019226 

0.41374813E-02 

407.060001 

416.131741 

-9.07173920 

-0.21S00162E-01 

414.990002 

428.071049 

-13.0810471 

-0.30558121E-01 

455.990002 

465.049557 

-9.05955505 

-0.19480838E-01 

489 . 7 7 OOOC 

480.385475 

9.38452530 

0. 19535406E-01 

598.510002 

600.017441 

-1.50743866 

-0.25123247E-02 

624.529999 

645. 126C60 

-20.5960617 

— 0.31925639E— 01 

1837. 89999 

1884.78838 

-46. 8883820 

-0.24877266E-01 

2386. 39999 

2349.77625 

36.6237488 

0.15586058E-01 

2554.00000 

2574. 14865 

-20. 1486511 

— 0.782 73067E-02 

2821.00000 

2775.28751 

45.7124939 

0.16471264E-01 

4349.00000 

4517.62238 

-168.622375 

-0. 37 3 2 5469 E -01 

5426.00000 

5254.61517 

171.384827 

0. 32616057E-01 

7009.00000 

6978.67584 

30.3241577 

0 • 43 4 52 5 95 E-02 

7587. 00000 

7467.71558 

119.284424 

0. 15973348E-01 

8788.00000 

8833.32520 

-45.3251953 

-0.51311589E-02 

9427. 50000 

9401.30420 

26. 1958008 

0.27864007E-02 

10360.0000 

10462. 5498 

-102.549805 

-0. 9801 607 2E— 02 

11060.0000 

11064.2125 

-4.21252441 

-0.38073423E— 03 

11720 .0000 

11846.6007 

-126.600708 

-0.10686670E— 01 

14020.0000 

13895.7180 

124.281982 

0. 8943905 OE -02 

14970.0000 

14985.5291 

-15. 5290527 

-0. 10362699E-02 


EQUATION FOR THE ABOVE IS 


(AO- A3) ARE 

0.68154073 -0 . 72963022E-05 -0 . 8688173 IE-1 1 

= 3157.9090 STANDARD DEVIATION' 56.195276 

.128715 
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FITLOS Fit 


SAMPLE PROBLEM 1 


SPLINE JOINTS CHOSEN BY PROGRAMMER 

DEGREE OF POLYNOMIAL = 2 NUMBER OF SEGMENTS = 3 

EQUATION FITTED IS Y * AO + A1 X + A2 X**2 

SEGMENT COEFFICIENTS IN ASCENDING ORDER - 


AO 

056095 51 19 1801 6D 00 
—1 •016499811 302 85 6 D 01 
-5.5008875652682030 01 


A 1 

5 .379181567783746 D- 01 
6.8012909302783340-01 
6 .92941 62 94677520 D- 01 


A2 

3.485214235920547D-04 
-7. 00591703161 387 ID- 06 
-7. 92 1098205894 84 OD- 06 


SPLINE JOINTS ARE 


8.8600000 200.00000 


7000.0000 47300.000 


X 

8. 8600000E+00 
1.09800006+01 
1 . 290Q000E+0 1 
1 • 68100006+01 
6. 0804 99 9E +01 
1 • 2086000E+02 
1 • 80880006+02 
%4082000E + 02 
2.4809000E+02 
2. 5589000E+02 
2.6594000E+02 
2.7578000E+02 
2 • 7930000E+02 
3.00850D0E+02 
3. 60890006+02 
3 . 747400 OE + 02 
4.0699000E+02 
4. 1579000E + 02 
4 . 455 3000E + 02 
4. 8048000E+02 
5.065 700 0E+02 
5.13070006+02 
5. 7347000E+02 
6.0079000E+02 
6. 26500006+02 
6. 442600 0E+02 
6.5930999E+02 
7. 22160006+02 
9. 0080000E+02 
9. 6833999E+02 
2.86540006+03 
3.59870006+03 
3.957 1000E + 0 3 
4,281 OOOOE+03 
7 • 200000 0E+03 
8. 5035999E+03 
1. 1750000E+04 
1.27300006+04 
1.5640000E+04 
1 • 6940000E+04 
1 • 95 5000 0E + 04 
2. 1160000E+04 
2.3440000E+04 
3.12300006+04 
4.73000006+04 


Y 

5. 5560000E+00 
7.6798C00E+00 
1.31160006+01 
1.1967000E+01 
4.30430006+01 
8.60379996+01 
1.16490006+02 
1.61030006+02 
1 .5800 000 E+02 
1. 5989000 E+02 
1.6800000 E+02 
1.74920006+02 
1 . 8184000 E+02 
2. 1051000 E+02 
2.4h47000E+02 
2. 380 3000E+02 
2 • 605 1 OCOE+02 
2. 7591000 E+02 
2.87060006+02 
3.20670006+02 
3.24270006+02 
3.2702000E+02 
3.72890006+02 
4. 0049000 E+02 
4. 0706000 E+02 
4. 1499000E+02 
4. 5599C00 E+02 
4.8 97 70 00 E+02 
5 • 985 1 COO E+02 
6. 2452999 E+02 
1.8379000E+03 
2.3864000E+03 
2.554C000E+03 
2.82100C0E+03 
4.3 49 OOOOE+03 
5.4260000E+03 
7.0089999 E+03 
7. 5 87 0000 F+ 03 
8 • 7880000 E + 03 
9. 42 74999 E+03 
1.0360000E+04 
1. 1060000 E+04 
1.1 72 OCOO E+04 
1.4020000E+04 
1.49700006+04 


Y* 

8. 84940918 1101 14 8D 00 
1 . 00 0445 47 65 98 7 3 5D 01 
1.105 32 3717141 96 5D 01 
1. 3 196 9 835 2055 7 19D 01 
3. 805277 90 80 25 758D 01 
7. 415978 48 14 54 96 2D 01 
1.12757 5067591 07 OD 02 
1.532173 86 8026619D 02 
1.58 1370239574 99 ID 02 
1 .634144897988261 D 02 
1 • 7021 3044896057 ID 02 
1. 768681700603464D 02 
1. 792485 36045032 3D 02 
1.938177278807263D 02 
2. 3437432806 94 36 9D 02 
2. 4372273793 86 33 OD 02 
2 • 6 548 02 7646 74 65 9D 02 
2 .71414685 881397 ID 02 
2. 9146226250 18356D 02 
3. 150060 349719130D 02 
3. 325701857941 790D 02 
3.3694459206281390 02 
3. 775646 1287 83 58 5D 02 
3. 959209841890134D 02 
4.1318604047068730 02 
4.2510702017359260 02 
4. 62029941 330561 3D 02 
4. 773433388 812 723D 02 
5. 96 8 10405 12 35 94 5D 02 
6.418618806031 82 9D 02 
1. 88 1 1 546991 33 1980 03 
2. 346684455466761 D 03 
2. 57 1470 70 41 79 836D 03 
2. 773070480931 50 3D 03 
4. 523541125521543D 03 
5. 264705 331817 194D 03 
6. 993448649542047D 03 
7.482501132021744D 03 
8.845021 5455 1 8303 D 03 
9 .41035547041391 ID 03 
1 • 046453644340334D 04 
1. 106099913502764D 04 
1. 183542541585308D 04 
1.3860008950 24 91 2D 04 
1. 4999336393 10552D 04 


DEV 

3. 29340 918682319 4D 00 
2 .32465479 2022 66 OD 00 
—2 .06276288 48 47I34D 00 
1. 229983512927797D 00 
-4. 99022 114 0994 86 2D 00 
—1.187 82 1529226191D 01 
-3.73249301201 1147D 00 
-7.812613883983602D 00 
1.37023957 499081 9D-01 
3. 52449040917763 1 D 00 
2 .21304489 60 57056D 00 
1.948169984052413D 00 
-2.591464107555623D 00 
-1.66922723 48 15548D 01 
-1.009567124391762D 01 
5.692737251987523D 00 
4.970278145932753D 00 
-4 .49531396 60 15033D 00 
4. 40226 112 854458 ID 00 
-5. 663 96 31 97032 302D 00 
8.3001853364152960 00 
9. 92459160 5050239D 00 
4.6746134887101220 00 
-4. 56901 748945341 5D 00 
6.1260390973963010 00 
1 .0117018495125820 01 
6.03993965209446 I'D 00 
-1.242666157649131D 01 
-1. 69959701 2635934D 00 
1.733188182388608D 01 
4.325 47 0523671356D 01 
-3.971553842972344D 01 
1.7470704179836500 01 
-4.7929519068497030 01 
1.745411255215436D 02 
-1.61294668 18280520 02 
-1.555135045795305D 01 
-1. 04498867978255 5D 02 
5.702154551830335D 01 
-1. 71445295 8608854D 01 
1.045364434033454D 02 
9.991350276395679D-01 
1. 154254158530803D 02 
-1.599910497508818D 02 
2.933639310551916D 01 


R-ERR 

0.5927662329400700 00 
0. 302697309813208D 00 
-0.1572707285756340 OC 
0. 102781274516891D 00 
-0. 115935718127078D 00 
-0. 138057780021801D 00 
-0.3204 13170172961D-01 
-0.4851651 152 37531 D-01 
0.8672 40237335961 D-03 
0.2204321985509880-01 
0. 1317 288628605390-01 
0. 1113749132862270-01 
-0. 1425134 24 185055D -01 
-0 • 79 2 9 44 3 888 65442 D-0 1 
-0.4129615606116680-01 
0.2391604938690780-0 1 
0. 1907903027890010-01 
-0. 1629268228226990-01 
0.1533568281015890-01 
-0. 1766 29033878 2 56D-01 
0.2559652550250760-01 
0.3034857681841400-01 
0. 1253617285623530-01 
-0. 1140856817974110-01 
0. 1504947446747160-01 
0.2437894516544150-01 
0. 1324 57 7212189270 -Ot 
-0. 2537244332008230-01 
-0.2839 71 3633 1 3838 D-02 
0.2775 188038647120-01 
0.2353485248400830-01 
-0. 1664244 82615240D-01 
0. 6840526303773100-02 
-0. 1699C2584432815D-01 
0.4013362279180120-01 
-0. 297262565762634D-01 
-0. 2218768791261670-02 
-0. 1377341083145580-01 
0.6488 569 13043962 0-0 2 
-0. 181856585373519D-02 
0. 1009 03902898982 0-0 i 
0.903 3770 593486 15 D-0 4 
0. 9848 584970399 34D-0 2 
-0. 1141 162979678190-01 
0. 195 9 67 8 898 1642 7 D-0 V 


CORRELATION OF FITTED DATA TO ORIGINAL DATA 


VARIANCE 
STANDARD DEVIATION 


3. 160794765434042D 03 
5. 62209458351135 2D 01 


CORRELATION INDEX 
MAXIMUM CORRELATION 


0.911037 79 41 55 017D OC 
0.911111111111111D CC 


NO REFIT CHECK MADE 


76 













APPENDIX G 


COMPUTER INPUT AND OUTPUT SHEETS FOR SAMPLE PROBLEM 2 


INPUT FOR SAMPLE PROBLEM 2 


Card column 

1-6 

7-12 

13-18 

19-24 

25-30 

31-36 


SAMPLE 

PROBL 

EM 2 



3 51 

0 0 

0 T 

F F 

F .0 

1 

(2E18 

. 7,F6. 

0) 




0 

. 0 


-1 

. 0 


6 

- 28318 

E-0 2 

-9 

.96054 

8 E-01 

1 

. 25663 

6 E— 0 1 

-9 

. 84250 

2 E-01 

1 

. 88495 

4 E-01 

— 9 

. 64679 

5 E-01 

2 

.51327 

2 E-01 

-9 

. 37497 

5 E-01 

3 

. 14159 

E-0 1 

“9 

.02919 

6 E-01 

3 

. 76990 

8 E-01 

-8 

.61220 

5 E-01 

4 

. 39822 

6 E-01 

-8 

. 12732 

8 E-01 

5 

. 02654 

4 E-0 1 

-7 

. 57844 

6 E-01 

5 

. 65486 

2 E-01 

-6 

. 96997 

|6 E-01 

6 

. 28318 

E-01 

-6 

. 30684 

2 E-01 

6 

.91149 

a E-01 

-5 

.59444 

9 E-01 

7 

.53981 

6 E-01 

-4 

. 83864 

4 E-01 

8 

. 16813 

4 E-01 

-4 

.04569 

1 E-01 

8 

. 79645 

1 E-01 

-3 

. 22222 

2 E-01 

9 

.42477 

E-01 

-2 

. 37520 

5 E-01 

1 

.00530 

9 

-1 

. 51190 

2 E-01 

1 

. 06814 

1 

-6 

. 39817 

5 E-02 

1 

. 13097 

2 

2 

. 33339 

4 E-02 

1 

. 19380 

4 

1 

.09970 

6 E-01 

1 

. 25663 

6 

1 

.95131 

4 E-01 

1 

. 31946 

8 

2 

. 78013 

9 E-01 

1 

.3823 


3 

.57815 

E-01 

1 

. 44513 

1 

4 

. 33736 

E-01 

1 

. 50796 

3 

5 

. 04987 

4 E-01 

l 

. 57079 

5 

5 

. 70795 

E-01 

1 

. 63362 

7 

6 

. 30403 

3 E-01 

1 

. 69645 

9 

6 

.83081 

8 E-01 

1 

. 75929 


7 

.2813 

E-01 

1 

. 82212 

2 

7 

. 64877 

6 E-01 

1 

. 88495 

4 

7 

. 92698 

7 E-01 

1 

. 94778 

6 

8 

. 11006 

6 E-01 

2 

.01061 

8 

8 

. 19262 

6 E-01 

2 

. 07344 

9 

8 

. 16979 

3 E-01 

2 

. 13628 

i 

8 

.03724 

E-01 

2 

. 19911 

3 

7 

.79122 

2 E-01 

2 

. 26194 

5 

7 

.42861 

2 E-01 

2 

. 32477 

7 

6 

. 94692 

3 E-01 

2 

. 3876C 

8 

6 

. 34434 

E-01 

2 

. 45044 


5 

.61973 

3 E-01 

2 

.51327 

2 

4 

. 77268 

5 E-01 

2 

. 57610 

4 

3 

. 8035C 

2 E-01 

2 

. 63893 

6 

2 

.71322 

1 E-01 

2 

. 70176 

7 

1 

. 50362 

2 E-01 

2 

. 76459 

9 

1 

. 77228 

9 E-02 

2 

. B2743 

i 

-1 

. 26269 

3 E-01 

2 

. 89026 

3 

-2 

.81214 

E-01 

2 

. 95309 

5 

-4 

. 46638 

E-01 

3 

.01592 

6 

-6 

.21996 

5 E— 0 1 

3 

. 07875 

8 

-8 

.C6675 

E-01 

3 

. 14159 


-i 

. 0 


2 51 

0 0 

0 F 

F F 

O 

* 

1 




FITLOS OUTPUT FOR SAMPLE PROBLEM 2 


sample problem 

DATA DIVIDED AS EVENLY AS 
DEGREE OF POLYNOMIAL = 
EQUATION FITTED IS 


POSSIBLE AMONG THE MAXIMUM NUMBER OF 
3 NUMBER OF SEGMENTS = 10 

Y = AO + A1 X + A2 X**2 + A3 X**3 


SUBSETS 


SEGMENT COEFFICIENTS IN ASCENDING ORDER - 


AO 

—9.9999926188 26976D— 0 1 
-9. 942762 028331 344D-01 
-9. 563760988676222D-01 
-8. 876375172 1469830-01 
-8. 8260209746658800-01 
-1.2012250572443010 00 
-2. 26 6 609 5 048 18 91 6D 00 
-4.5532967150211330 00 
-8.5195615738630290 00 
-1.3815877854824070 01 


A 1 

-1.2196927085028620-03 
-5.5870927477371880-02 
-2.36830747977 3648 D-01 
-4. 556326321 1806770-01 
-4 .67653831467032 4D-01 
1 .408 71 778130531 3 D-01 
1 .836485221982002D 00 
4 .9559525400400 1 60 00 
9 .690 3365254402 1 6D 00 
1 .5309907317161560 01 


A2 

1.022398505 109 77 70 00 
1.196358936082 10 6D 00 
1.4843656807206570 00 
1.7165218839 94 57 40 00 
l. 726088050 75 2 87 80 00 
1.338688 30 2785158D 00 
4. 39 1 366 541 3 8565 ID- 01 
-9. 79374960064888 OD-Ol 
-2.8631 28 043711185D 00 
—4. 850646 227598 1900 00 


A3 

-1. 1231313538701 22D-01 
-2.96891042 17686 69D-01 
-4.4968349530245 18D-01 
-5.3179202490355 240-01 
-5. 3432953 19471508D-01 
-4. 5212068222463 13D-01 
-2.9304489120841 03 D-01 
-7. 8032205 8 1263 304D-02 
1.71808516 14102 72D-01 
4. 061 2229146063 33D-01 


SPLINE JOINTS ARF 


0 

2.1991130 


0.3141590 

2.5132720 


0.6283180 

2.8274310 


0. 9424770 
3. 1415900 


1.2566360 


1.570795C 


X 

0. 

6 . 283 180 0E-02 
1. 2566360E-01 
1. 8849540E-01 
2 • 513272 OE-0 1 
3. 1415900F-0 1 
3. 7699080E-01 
4 . 3982260E-0 I 
5. 0265440E-0 1 
5.6548620E-01 
6.2831799E— 01 
6.9114980F-01 
7.5398160E-01 
8.168 1340E— 0 1 
8 • 796451 OE-O 1 
9.4247700E— 01 
1 . 0053090E+00 
1 • 068 141 OE+OO 

1. 1309720E + 00 
1. 1938040E+00 
1. 25 66 36 OE+OO 
1. 3 19468 OE+OO 
1 . 382 30006+00 
1 • 445 1 31 OE+OO 
1 . 5079630E+00 
1 • 5707950E+00 
1 • 6336270E+00 
1 .6964590E+00 
1 ■ 7592900E+00 
1.822 1220E+00 
1.8849540E+00 
l • 9477860E+00 
2.0106180E+00 
2 .073449 OE+OO 
2 • 136281 OE+OO 

2. 1991130E+00 
2.261 9450E+00 
2.3247 77 OE+OO 
2 • 387608 OE+OO 
2 .450440 OE+OO 
2 . 51 3272 OE + OO 
2 • 576 1040E+00 
2 . 6389360E+00 
2. 701 7670E+00 
2.7645990E+00 
2. 82 743 1 OE+OO 
2 . 8902630E+00 
2 . 9530950E+00 
3.0159260E+00 
3 • 078 7580E + 00 
3.141 5900E+00 


Y 


-1 

. 000 OC 00 E+00 

-9 

• 960 5479 E- 

-01 

-9 

. 842 5020F- 

-01 

-9 

. 646 7949 E- 

-01 

-9 

. 3749750 E- 

-01 

-9 

.0291959E- 

-01 

-8 

.612 2049E- 

-01 

-8 

. 1273279E- 

•01 

-7. 

.5784460 F- 

-01 

-6, 

. 9699759 E- 

-01 

-6, 

. 306 8420 E- 

•01 

-5, 

. 5944489E- 

-01 

-4, 

, 8 38 6440 6- 

-01 

-4. 

• 045691QE- 

-01 

-3. 

.22222206- 

01 

-2. 

, 3752050E- 

■01 

-1, 

.5119020 E- 

•01 

-6, 

. 398 1750E- 

■02 

2. 

.3333940E- 

■02 

1 ■ 

► 099 7060 F- 

01 

1. 

.9513140 E- 

■01 

2, 

.7801390E- 

■01 

3. 

,5781 500 E- 

01 

4, 

, 3373600 E- 

■01 

5, 

,049 8 739E- 

01 

5. 

, 707 9 500 E- 

01 

6, 

, 3040330E- 

01 

6. 

.8308 179 E- 

01 

7. 

,281 2999E- 

■01 

7, 

, 648 7 760E- 

01 

7, 

,92698706- 

01 

8. 

,11006596- 

01 

8, 

, 192 6260 E- 

01 

8. 

, 1697929E- 

01 

8, 

.0372399E- 

01 

7. 

. 791 2 220 E- 

01 

7 . 

, 42861 19E- 

01 

6. 

, 9469230 E- 

01 

6. 

.344 3400 F- 

01 

5. 

.6197330 E- 

01 

4. 

,77268496- 

01 

3 « 

, 8035C20 E- 

01 

2. 

, 7 13 2210 F- 

01 

1. 

,50362206- 

01 

1. 

, 7722 890 E- 

02 

-1. 

,2626930 6- 

01 

-2. 

■8121400E- 

01 

-4. 

4663799E- 

01 

-6. 

,2199650 E- 

01 

-8. 

,066 7499E- 

01 

-1. 

•OOOOOOOE+ 

00 


Y* 

-9. 99999261 8826976D-01 
-9. 960674958476 492D-01 
-9. 84230363673541 5D-01 
-9. 6465 5021 41 6831 9D-01 
-9.375086223605473D-01 
-9. 02958326 21 25 283D-01 
-8. 61 2170680242 639D-01 
-8. 1268092733 52 72 2 D-01 
-7. 5779176293 12658D-01 
-6. 9699144732 18704D-01 
-6. 307218 364953882D-01 
-5. 59462683 7207597D-01 
-4. 8 38453690 1 71 605D-01 
-4. 045391 37 1470 360D-01 
-3.222133985974475D-01 
-2. 3753701006096250-01 
-1.51 1997302 793477D-01 
-6. 3 972 79 22 7 545 349D -02 
2.3350947147213090-02 
1. 09982801 3747202D-01 
1.95 1299093027 130D-01 
2. 780001 7472596590-01 
3. 577989657881 044D-01 
4. 337298688772676 D-01 
5 .0500005876655920-01 
5.708130472710383D-01 
6. 303940 399477 847D-01 
6. 83049625323915 ID-01 
7. 28 10631 56077 809D-0 l 
7.6489265515878380-01 
7.9273497492909770-01 
8. 1099977515802130-01 
8. 1921155793454400-01 
8. 1693425828156660-01 
8. 0373165758366840-01 
7.791675770659641D-01 
7.4285894970995780-01 
6.946368094941 11 8D-01 
6. 343858 40 2942 706D-01 
5. 619879640946 133D-01 
4. 773279 190 17 2 62 ID-01 
3. 80 35 17446484 16 8D-01 
2. 7125271 751 28763D-01 
1 • 5028888 84278804D -01 
1 .7712008085193350-02 
-1.262203631222496D-01 
-2. 81 194190911961 ID -01 
-4. 46663576 1433 189D-01 
-6. 220210562781 592D-01 
-8. 06667692 3403300D-01 
-9. 999963788338 42 OD-Ol 


DEV 

7.3 8 117 302 3901 51 ID-07 
-1.269748 30 0997 13 ID-05 
1. 98391 0145997303D-05 
2. 44805339089 1047 D-05 
-1 .11257 55627949820-05 
-3. 872828877271584 D-05 
3.4333390135365960-06 
5.187350605917 46 2D-05 
5.283416994572576 D-05 
6. 15049173 580 1289D-06 
-3. 763954638308 364D-05 
-1. 77868439 4776464D-05 
1.90315 18661 44981 D-05 
2.99634415241 1679D-05 
8. 8 00 69 3543309 36 8 D- 06 
-1. 6510906037 15932D-05 
—9. 52945 878 25200950-06 
8. 957 77 32 4 350088 8D-06 
1. 70071 907 1195955D-05 
1. 22010304794972 8D -05 
-1 .4912973627101420-06 
-1.37251 9640496 78 3D-05 
- 1.60353 718 0194792D-05 
-6. 1313327791 5 78940-06 
1.266 24667 201 591 4D-05 
1. 80476715 8153364D-05 
-9. 262662 13 574500 3D-06 
-3. 2172931528 46 246D -05 
-2. 36822416 834 758 8D-05 
1. 50527007836 1320D-05 
3.627122283322848D-05 
-6. 8230 1 655063 8883D-06 
-5.1 0435006 4044638D-05 
-4. 503822 38 8457479D-05 
7.65921 5881 972 1590-06 
4.537347747834986D-05 
-2.24668 877 72 40990 D-06 
-5. 54892758 3909347D-05 
-4.815936336033388D-05 
1.46605381 7100743D-0 5 
5.942062634911 59 3D-05 
1.54605 l 70 0884732D-06 
-6. 938 38 4172272549 D-05 
-7. 33107444 97 23 47 4D-05 
- 1.088 19286 355 1996D-05 
4. 89364146938342 5D-05 
1.98078 825 945 9 132D-05 
-2. 5576876932 83128D-05 
-2. 45566 801 3290 33 5D-05 
7. 309638544228392D-06 
3.621 16615 804097 8D-06 


R-ERR 

-0. 7381 173023901 51 D-06 
0. 127477755549379D-04 
-0. 201 5 656 324381 60 0-0 4 
-0. 2537 685 71420941 D-04 
0. 118675043594686D-04 
0.4289 228 948 15560 D-04 
-0.398659693783620D-05 
-0.6382602745388860-04 
-0.6971 636420952 10D-04 
-0.8824 265 32759 167D-05 
0.5968049709374130-04 
0.3179 373705446680-04 
-0.3933 233 906 10565 D-04 
-0.7406 260 507914 84 D-04 
-0. 2731 249 91470922D-04 
0.6951 360 44926796 D-04 
0.6302 960 595858 85D-0 4 
-0. 1400051301610700-03 
0.7288 6065292291 7 D-03 
0. 1109481119617830-03 
-0.7642528870925160-05 
-9.4936 87416665148D-0 4 
-0.4481 469 96351 943 D-04 
-0. 141 3 60 9 379020 57D-0 4 
0.2507 481 733 78979D-04 
0.3161 84822820771 D-04 
-0. 146932322489378D-04 
-0. 47099676 21 83248 D-04 
-0.325 2 47438691 188D-0 4 
0. 19679881 7683248 D-04 
0. 457566319506293D-04 
-0. 8413022244203990-05 
-0.6230419959488980-04 
-0. 551277420091566D-04 
0 *9529659307830300-05 
0.582 3 666 33493031 D-04 
-0.30 24372235530 08 D-05 
-0. 7987604862950820-04 
-0.7590 9 17792287 75D-04 
0.2608 760 60806241D-04 
0. 1245014631 16568D-0 3 
0. 406481 1 078 28737D-0 5 
-0.2557 24 9 90 24 70 93D -03 
-0.4875610020387510-03 
-0.6140041 848 15736D-0 3 
-0 • 387555921 1 40 725D -03 
-0.7043 704 324642 18 D-04 
0.5726 53401 072948 D-04 
0. 394804153219085D-04 
-0.9061441753242550-05 
-0.3621 166 15804098D -05 


CORRELATION OF 


FITTED DATA TO ORIGINAL DATA 


VARIANCE 
STANDARD DEVIATION 


1.9405527708483540-09 
4. 40 51705572 20 11 7D-05 


CORRELATION INDEX 
MAXIMUM CORRELATION 


0.470588 232 238856D OC 
0.470588235 29 4118D CC 


REFIT CHECK WAS MADE 

DUPLICATION OCCUREO IN FIRST SET OF COEFFICIENTS - CURVE WAS REFIT IN NEW SEGMENTS 


1.8849540 


82 



SAMPLE PROBLEM 2 

DUPLICATION OCCURED IN FIRST SET OF COEFFICIENTS - CURVE HAS REFIT IN NEW SEGMENTS 


DEGREE OF POLYNOMIAL = 


NUMBER OF SEGMENTS = 


EQUATION FITTED IS Y 

SEGMENT COEFFICIENTS IN ASCEND 
AO 

-1.0000533713767180 00 -6 

-9.2222092697920740-01 -3 

-9.59871B8695653460-01 -2 

—4. 228361 103683710D 00 4 

—4. 694865703582 764 D 00 5 

— 1.2000 20 6 30 65767 3D 01 1. 


» AO + A1 X + A2 X**2 + A3 X**3 
ING ORDER - 


A 1 

.6644727526181670-03 1, 
.7827301246579740-01 1 
.72525B497230243D-01 1 
.6043168194592000 00 -8, 
.1754450 79 80346 7D 00 -1, 
.3480314329266550 01 -4 


A2 

0764140970097740 00 
6678479594265810 00 
5688468323423880 00 
5669728182256220-01 
089768171310 42 5D 00 
2368 20 9129 19 93 9D 00 


A3 

-2. 1091791470 87 161D -01 
-5.246835903 1990 99D-01 
-4. 937 8844099192 07D-01 
—9. 166595758 8702 44D-02 
-5. 9961467981338 50 D-02 
3. 3 75 5400 299 56 028D -01 


SPLINE JOINTS ARE - 


X 

0. 

6.2831 80 OE-02 
1 . 2566360E-0 1 
1 . 8849540E-01 
2.5132720E-01 
3.14 1 5900E-0 1 
3. 76990B0E-01 
4. 3982260E— 01 
5. 026544 0E-01 
5.654 862 OE-0 1 
6. 2831799E— 01 
6.91 14980E— 01 
7. 5398160E-0 1 
8. 168 I340E-0I 
8. 79645106-01 
9.4247700E-01 
l • 005 3090E + 00 
1 . 068 141 OE + OO 
1. 1309720E+00 
1. L938040E+00 
1 • 2566360E+00 
l • 319463 OE+OO 
1. 3823000E+00 
1.445131DE+00 
1. 5079630E+00 
1.5707950E+00 
1 . 63362706+00 
1.6964590E+00 
1 .75929006+00 
1 .8221220E+00 
1.88 49 540 E + 00 
1 • 9477860E+00 
2 • C 106 180E+00 
2.0734490E+00 
2. 136281 OE+OO 
2.19911306+00 
2.261 945 OE + OO 
2. 32477706+00 
2.3876080E+00 
2 . 45044006+00 
2.5 1 32720E +00 
2. 5761040E+00 
2.6389360E+00 
2.701 7670E+00 
2 . 764 599 OE + 00 
2 • 027431OE+OO 
2 . 890263 OE + OO 
2.9530950E+00 
3.01592606+00 
3.0787580E+00 
3 ■ 14 1 5900E+00 


0 0 


Y 


-1 . 0000 COO E+ 00 

-1, 

-9.960 5479 E- 01 

-9, 

-9.8425020E-01 

-9, 

-9. 646 7 949 E- 01 

-9, 

-9. 3749750 E- 01 

-9, 

-9 .029 1 959 E-0 1 

-9, 

-8. 6122049 E— 01 

-8, 

-8.1273279E-01 

-8, 

-7. 5 784 A 60 6-0 1 

-7. 

-6 .9699 7 59 E-0 1 

-6, 

-6 . 306 8 420 £-0 1 

-6, 

-5. 5944489 E- 01 

-5, 

-4. 8386440 E-0 1 

-4, 

-4 .0456910 E- 01 

-4, 

-3. 22222206-01 

-3. 

-2.37520506-01 

-2, 

-1. 5119020 E-01 

-1, 

-6. 3 98 1750 E- 02 

-6. 

2. 3333940 E-02 

2. 

1.0997060 E-01 

1 

1.9513140 E-01 

I 

2.7801 390E-0 1 

2, 

3. 5781500 E-01 

3. 

4.337 3600E-0 i 

4, 

5. 04987 39 E-01 

5, 

5.7079500E-01 

5. 

6.3040330E-01 

6, 

6 » 8308 179E-0 1 

6, 

7.2812999E-01 

7, 

7 . 648 7 760 E-0 1 

7. 

7.9269870E-01 

7. 

8.U00659E-01 

8. 

8 . 192 6 260 E-0 1 

8, 

8.169 7929E-0 1 

8, 

8 .037 2 399E-01 

8, 

7. 7912220 E-01 

7, 

7.4286 119E-01 

7, 

6. 9469230 t-01 

6 1 

6.344 3400 E-01 

6 

5. 6197330 E-Oi 

5 

4. 7726849 E-01 

4, 

3.8035020E-01 

3 

2.7132210E-01 

2 

1.5036220F-01 

1 , 

1.7722090 E-02 

1 

-1 .2 62 6930 E-01 

-1 

-2. 81214 00 E-01 

-2 

-4.4663799E-01 

-4 

-6. 2199650 E-01 

-6 

-8. 066 7499 E-01 

-8 

-l.GOOOOOOE+OO 

-9 


6203180 

1.0681410 

Y* 



00005037 1376718D 00 -5. 
96 2719248977712 D-01 -2. 
8430837622605160— 01 -5. 
644736346767 1 79D-01 2. 
3708 1606 758655BD-01 4. 
0244620545221630-01 4, 
6088133454 68 538D -01 3, 
1270090449806070-01 3, 
58 21882050 85 834 D-01 -3. 
9774900142454650-01 -7. 
3160534971065350-01 -9. 
6017959462390150-01 -7. 
8417480885714940-01 -3. 
0437185976867630-01 l, 
2155178111428740-01 6. 
3649504982647670-01 1. 
4998268161816150-01 1, 
2795698742343480-02 1. 
4291290648083220-02 9, 
105384892889192D-01 5. 
9520959951546090-01 7. 
7756972947798850-01 -4. 
5688392891639910-01 -9, 
3241614034791210-01 -l, 
0343386285108560-01 -1, 
6920095319309130-01 -1. 
2898250032391780-01 -1. 
82043593 1935 54 5D-01 -1 . 
2764865921640740-01 -4. 
650642463928 135D -01 1. 
9355466347578140-01 8. 
1 2 384 98 5 7 56 554 4D -01 l, 
2082031002926390-01 1. 
182256 18 7695032D-01 1. 
0436461848796090-01 6, 
7910086014383140-01 -2. 
4229791234001220-01 -5, 
9381934757122380-01 -8, 
3352980660373690-01 -9. 
6129091917981100-01 -6 
7697677278847730-01 -2. 
804866687487092D-01 1 
7173312200208110-01 4 
507260357413340D-01 3 
7868523246207250-02 1 
2633684850429280-01 -6. 
813876939819000D-01 -1 
467816276583525D-01 -1 
220 13446 75 5247 8D- 01 -1, 
0658625800850810-01 8, 
9999500873941290-01 4 


2.0136180 2.4504400 2.638936C 


DEV 

03713767 1770 60 8D-05 
17126533131 994 9 D- 04 
8173451350080870-05 
05867 2740229417D-04 
158 898462636 174D-04 
73392471 53932450-04 
3916681642365280-04 
1 8963 432 707 208 8D-05 
7422340737192080-04 
5140361094028220-04 
211527616483917D-04 
3469774708961880-04 
10408321 327448 9D-04 
9724081988375940-04 
704181 76 7034356D-04 
0254493284486690-03 
2075192024036950-03 
186051 33635351 9D-03 
5735069158209280-04 
6788894467851850-04 
8198915385185860-05 
4417044438238570-04 
3107224350730270-04 
3198598621346760-03 
55353344875341 9D-03 
59404640636540 ID-03 
420 80228 600260 90-03 
03820506 18891 19D-03 
8133863305699350-04 
8664393481337970-04 
5595976951691190-04 
3783875819823880-03 
5577085940794430-03 
2463222640519690-03 
40620 120 17451760-04 
1343444654320360-05 
6328405872285710-04 
72951 19872710100-04 
,0419305389399800-04 
,8238437663126210-04 
917256024357329D-04 
, 36470151 99 32510D-04 
,1102064748210760-04 
,6383656895644380-04 
4563323237837270-04 
, 7548967349395640-05 
,7369518734433600-04 
,43628391 9664660D-04 
.694715 72 2155848D-05 
.8743970366067510-05 
.9912605870616740-06 


R-ERR 

0.5037137671770610-04 
0. 2179865339622690-03 
D. 5910433229890760-04 
-0.2134048392306920-03 
-0.4436170205997700-03 
-0. 5242908367786900-03 
-0.3938 211130445280-03 
-0. 3924579300564980-04 
0. 4937996639460670-03 
0. 1078057676665370-02 
0. 1460561032137090-02 
0.1313262040982380-02 
0. 6415 192 37587450 0-03 
-0.4875 330 804967970-03 
-0.2080608282665210-0? 
-0.4317308746390800-02 
-0. 7986755727884750-02 
-0. 1853733768536020-01 
0.4102824869553860-01 
0.516400695186584D-02 
0.400 7 500 337962 290-03 
-0.1597655529117110-02 
-0.2602105111549560-02 
-0. 3043002797774460-02 
-0.3076380638678360-02 
-0. 2792 677 594379 75 D-02 
-0.22 5 3 79 8 925418650-02 
-0. 1519883950268680-02 
-0.6610613962872420-03 
0.2440180418587020-03 
0. 1079804679274570-02 
0. 1699600946632110-02 
0. 1901 354451369620-02 
0. 1525524905444960-02 
0. 7970648 151299420-0 i 
-0.2739421948959560-04 
-0.7582628645210740-03 
-0. 1256601232333770-02 
-0.1425 196 40243420 0-02 
-0. 1214264756551240-02 
-0.6112400114804790-03 
0. 3588013 165150200-03 
0. 1514880820360930-02 
0 . 2419734287999710-02 
0. 8217239528357820-02 
0. 5349595475472830-03 
0.6176619517125500-03 
0. 3215767404528480-03 
0. 2724638680846410-04 
-0.1100 120496462240-03 
-0.4991260587061670-05 


CORRELAT ION OF 


FITTED DATA TO ORIGINAL DATA 


variance 

STANDARD DEVIATION 


8. 1663887434215320-07 
9. 0368073870195080-04 


CORRELATION INDEX 
MAXIMUM CORRELATION 


0. 705881 863095631D OC 
0.705882 352941 1760 OC 


NO REFIT CHECK MADE 


3.1415900 
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SAMPLE PROBLEM 2 


DATA DIVIDED AS EVENLY AS POSSIBLE AMONG THE MAXIMUM NUMBER OF SUBSETS 
DEGREE OF POLYNOMIAL = 2 NUMBER OF SEGMENTS = 10 

EQUATION FITTED IS Y = AO + A1 X + A2 X*»2 

SEGMENT COEFFICIENTS IN ASCENDING ORDER - 
AO 

-9. 9 99730 3068974430-01 
-1.0185126112251060 00 
-1.16018*0887528850 00 
-1 .5752865214463010 00 
—2 .3800351341269560 00 
-3.5423*88166634340 00 
-4.8097794992208950 00 
-5.6720548603625500 00 
-5.4018257222269310 00 
—3. 22 17 16 79 8 6 673 30D 00 


A1 

l .505 1882719490100-03 
1.1953192395549190-01 
5.7048662598936060-01 
1.4513620881729570 00 
2.732 1603729542 40D 00 
4.2120653673700870 00 
5.5568522204412150 00 
6.3410550845583200 00 
6.1 260 1338 6 3 5431 60 00 
4 .583900472 55 6 707 D 00 


A2 

9. 8019617612206390-01 
7. 9235 06 325062 13 5D-01 
4.3349197790189460-01 
-3.3827362316472430-02 
-5.4344124063300110-01 
-1.0145100289810220 00 
-1.3712261210203000 00 
-1.5495259310300750 00 
-1. 506744 708122824D 00 
-1.2340390404242500 00 


SPLINE JOINTS ARE - 

0 0.3141590 0.6283180 0.9424770 1.2566360 1.570795C 

2.1991130 2.5132720 2.8274310 3.1415900 


1.8849540 


X Y Y* DEV 

0. -1 .00000006+00 -9.9997303068974430-01 2.6969310255697820-05 

6 • 283 1800E-02 -9.960 5479E-01 -9.9600880408714610-01 4.5994277493133230-05 

1.2566360E-01 -9.8425020E-01 -9.8430527165758000-01 -5.5068882578501730-05 

1. 8849540E-01 -9 . 646 7949F-01 -9.6486243409214330-01 -1.8293214140258090-04 

2. 5132720F-01 -9.3749750E-01 -9.3768028931754380-01 -1.8279271262444660-04 

3. 14159D0E-01 -9.029 1959E-01 -9.0275884170699390-01 1.6075621676170380-04 

3. 7699080E-0 1 -0 . 6 122049 E-0 1 -8. 6083966966767030-01 3. 80831 695 60709780-04 

4. 398 2260E-01 -0 . 1273 279 E-0 1 -8.1266435827911880-01 6.84425622 12611690-05 

5 • 0265440E-0 1 -7 . 57844 60E-01 -7.5823290412888980-01 -3.8830702767833620-04 

5.6548620E— 01 -6. 9699759E-01 -6.9754531745507570-01 -5 . 477 19641 469490 7D-04 

6.283 1799E-0 1 -6. 3068420E-01 -6. 306015921741 3290-01 8.2604774872208080-05 

6.91 14980E-0 1 -5 . 5944489E-01 -5 . 588184342 1 142880-01 6 . 26 462 665 38 3 1 3 1 8D-04 

7. 5398 160E-0 i -4 . 838 6440 E-01 -4 . 83612 57523299060-01 2.5182530283129890-04 

8. 1681340E-01 -4 .045 69 10 E-0 1 -4.0498399740250380-01 -4.14896813 94368820-04 

8 • 79 6451 OE-O 1 -3. 222 2220fc-0 1 -3.2293284849068940-01 -7. 1 0649 199698 543 1 D-04 

9 . 424 7700E-0 1 -2 . 375 2C50E-0 1 -2.3745872166868900-01 6.17774862 36369580-05 

1.0053090E+00 -1. 511 9020E-01 -1.5040664579161330-01 7 . 83555 028 95 1 88780-04 

1 .0681410E+00 -6 . 398 1 750 E-0 2 -6. 362 166 1292708960-02 3.60088 75598803060-04 

1. 1309720E+00 2 . 333 3940E-02 2 . 289485921 1 76881 D-02 -4. 390 83 744 732 3 1200-04 

1. 1938040E+00 1 . 099 7060 E-01 1 .09 145665 1 981 644D-0 1 -8. 24935 14607628880-04 

1 . 256636 OE + OO 1 .951 3 140E-01 1 . 95 129379 8063 12 70-01 -2 . 02079 376299 1 246D-06 

1. 3194680E+00 2.78013906-01 2.7883413814864350-01 8.2023822617260930-04 

1. 3823000E+00 3 . 578 1 500E-0 1 3. 5824801631203030-01 4. 330 15 152 12390560-04 

1.4451310E+00 4. 3373600 E-01 4. 333698750621 7420-01 -3.6612514787259800-04 

1.5079630F + 00 5 . 049 8 7 39 E-0 1 5 . 04 202 101 75542880-0 1 - 7. 85294544410231 ID-04 

1.570 7950E+00 5 . 7079 500 E- 0 1 5.7074346862346000-01 -5. 1530975996749720-05 

1. 6336270E+00 6 . 3040 330 F-0 1 6.3113426200474490-01 7.3095947482437750-04 

1. 6964590E+00 6 . 830 0 1 79 E-0 1 6. 03514768483921 10-01 4. 3297022847754720-04 

1.759290 OE+OO 7 . 28 1 2 999 E-0 1 7 . 278843464577922D-01 -2 .456 51 391 672 208 1 D-04 

1 . 822 1220E+00 7. 6487760E-01 7. 6424440606196730-01 -6.3319639603287570-04 

1 • 8849 540E+ 00 7 . 926 9 8 70 E-0 1 7.9259418440189890-01 -1 . 045 19304 36 56 1 34D-04 

1 ■ 947 7860E + 00 8 . 1 1006 59 E-0 1 8 . 1 1 525402958561 70-01 5.18804783 98969030-04 

2.01061806+00 8 . 1 92 62 60 E- 0 1 8 . 19629e040262473D-01 3.67202 59106282070-04 

2 • 0734490E+0C 8 . 1 69 7 9 29 E-0 1 8.1690751890833550-01 -7.1777597115740830-05 

2. 136281 OE + OO 8 . 037 2 399 E-0 1 8.03 35845080569520-01 -3.65547 56209117300-04 

2.1991 130E + 00 7 • 791 2220 fc-9 1 7.7098256917293910-01 - 1 . 396 344 1 5 546652 2D-04 

2.261 9450E + 00 7 . 428 6 1 19 E- 0 1 7.4307596720976580-01 2.14773 011030 71780-04 

2.3247770F+00 6 . 946 92 30 t~0 1 6.9493474275128090-01 2 . 42443 98 1 329 99830-04 

2 • 38 7608 0F + 00 6 . 344 34 00 E-0 1 6.3455996786358780-01 1.2596820595689450-04 

2 .4504 43 OE+OO 5 . 6 19 7 3 30 E-0 1 5. 6194969521 163700-01 -2.3638344805303720-05 

2.51 32720E+00 4. 772 6 849E-0 1 4.7710479981030750-01 - 1 . 63698 5 80605 4980D-04 

2 . 576 104 OE + OO 3 . 803 5 020 E-0 1 3.8019417600890650-01 -1.5 6022 58 7 809445 70-04 

2. 630936 OE + OO 2 . 7 1 3 2 2 10 E-0 1 2.7138666351090900-01 6. 45621 563 1001072D-05 

2. 7017670E+00 1 . 5036220 E-0 l 1.5068440810625190-01 3 . 222089338743572D-04 

2.7645990E+00 1 . 7722 890 E-02 1 . 8083467862419770-02 3. 60577848 59009280-04 

2. 827431 OE+OO - 1 . 2626 930 E-0 1 -1.2641430657717790-01 -1 .4500704023445720-04 

2 • 890 2630 F + OQ -2 . 8 12 1 400 E-0 1 -2.8173231168454380-01 -5.18312 889 98 815 95 0-04 

2 . 95 3095 OE + OO -4.4663 799E-0 1 -4.46 79394315 86 8910-01 - 1 . 55943 892 303 1 1200-04 

3.01 5 926 OE + 00 -6. 21 99650 E-0 1 -6.2159638886704680-01 4. 00 1 1073097945270-04 

3. 078 7560E + 00 -8 . 066 7499 E-01 -8.0614512107911460-01 5. 29880899759582 7D-04 

3. 141 590 OE + OO -1 .OOOOOOOE + OO -1.0004374799161070 00 -4 . 374799 16 1066534D-04 


R-ERR 

-0.2696931025569780-04 
-0. 46 1764529 1C67 160-04 
0. 5595 C08 507312490-04 
0. 189629966255800D-03 
0.1949794141172830-03 
-0. 1780404557962410-03 
-0. 4421999882774000-03 
-0. 8421 287062828120-04 
0.5123834479570450-03 
0.7858 271 580671410-03 
-0. 1309764463289500-03 
-0. 1119793332427300-02 
-0.5204460227957100-03 
0. 1025527686963990-02 
0. 2205463190500950-02 
-0.2600 932022900250-03 
-0.51825781 34688920-02 
-0.5627991 665C2892D-02 
-0.1881725698921150-01 
-0.7501 415319130720-02 
-0. 1035606651096040-04 
0.2950349699787100-02 
0. 1210164891690480-02 
-0.8441 198048935140-03 
-0.1555077513150360-02 
-0.9027 930 52 3727520 -04 
0. 1159510858839330-02 
0.6338482801308380-03 
-0.3373729861394810-03 
-0.8278401 589980470-03 
-0. 1318524981521140-03 
0. 639704763385927D-03 
0.4482 110991268920-0 3 
-0.8785730240994170-04 
-0.4548 17 2 79107669 0-0 3 
-0. 1792 201722 70182 0-03 
0. 2891 129 757105 230-0 3 
0. 3489 947 733 108 30 D-0 3 
0. 1985521047498600-03 
-0.4200 9726 54020 85 0-0 4 
-0. 3429905412936320-03 
-0.4102 077201092140-0 3 
0.2379 539 152456 750-0 3 
0.2142 885 21 74805 7 C-02 
0. 2034531 886703360-01 
0. 1148395063299070-02 
0.1843 126203567200-0 ' 
0.3491505258380470-03 
-0. 6432 684 600965 3 2 C-0 3 
-0 .6568703610000540-0 3 
0.4374799161066530-03 


CORRELATION OF FITTED DATA TO ORIGINAL DATA 

VARIANCE = 2.4801717257076730-07 CORRELATION INDEX = 0.6470587006306660 OC 

STANDARD DEVIATION = 4.9801322165876630-04 MAXIMUM CORRELATION = 0.6470588235294120 CC 

NO REFIT CHECK MADE 
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TABLE I. - COMPARISON OF FITLOS CURVE WITH f(x) 


0 

0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 

1. 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 

2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
2 . 
3. 
3. 
3. 


X 

y 

y 

y -y 


-1.0000000 

-1.0000504 

-0.5036592E— 04 

6283180F-01 

-0.9960548 

-0.9962719 

-0.2171248E— 03 

1256636 

-0.9842502 

-0.9843084 

-0.5816668E-04 

1884954 

-0.9646795 

-0.9644736 

0.2058744E— 03 

2513272 

-0.9374975 

-0.9370816 

0.4158914E-03 

3141590 

-0.9029196 

-0.9024462 

0.4733950E-03 

3769908 

-0.8612205 

-C. 8608813 

0 . 3 391 72 8E -03 

4398226 

-0.8 127328 

-0.8127009 

0. 3190339E-04 

5026544 

-0.7578446 

-0.7582188 

-0.3742203E-03 

5654862 

-0.6969976 

-0.6977490 

-0.7513985E— 03 

6283180 

-0.6306842 

-0.6316053 

-0.9211 525E-03 

6911498 

-0.5594449 

-0.5601796 

-0.7346943E-03 

7539816 

-0.4838644 

-0.4841748 

-0.3104061E-03 

8168134 

-0.4045691 

-0.4043719 

0.1972429E-03 

8796451 

-0.3222222 

-0.3215518 

0.670421 9 E— 0 3 

9424770 

-0.2375205 

-0.2364950 

0. 1025449E-02 

0053090 

-0.1511902 

-0.1499827 

0.120751 9 E— 0 2 

068 1410 

-0.6398175E-01 

— 0.6279570E-01 

0.118&051E-02 

1309720 

0.23333948-01 

0.242 91 29 E-0 1 

0.9 573 50 5 E-03 

1938040 

0. 1099706 

0.1105385 

0 • 5678888E-03 

2566360 

C. 1951314 

0.1952096 

0.7819757E-04 

3194680 

0.2780139 

0.2775697 

-0.4441738E-03 

3823000 

0.3578150 

0.3568839 

-0.9310730E-03 

4451310 

0 .4337360 

0.4324161 

-0.1319863E-02 

5079630 

0.5049874 

0. 5034339 

-0. 1 553535E— 02 

5707950 

0.5707950 

0.5692009 

-0.1594052E-02 

6336270 

0.6304033 

0.6289825 

-0.1423803E-02 

6964590 

0.6830818 

0.6820436 

-0.1038209E-02 

7592900 

0.7281300 

0. 7276487 

-0 .4813448E-03 

8221220 

0.7648776 

0.7650642 

0. 1866370E-03 

B849540 

0.7926987 

0.7935547 

0.8559 52 5E-03 

9477860 

0.8110066 

0.8123850 

0.1378387E-02 

0106180 

0.8 192626 

0.8208203 

0.155771 5E-02 

3734490 

0.8 169793 

0.8182256 

0.1 24631 8E-02 

1 36 28 1 C 

0 . 8 C 3 72 40 

0.8043646 

0.6406158E-03 

1991130 

0.7791222 

0. 7791009 

-0 .2134591E-04 

2619450 

0.7428612 

0.7422979 

-0.5632862E-03 

3247770 

0.6946923 

0.6938193 

-0.8729547E-03 

3876080 

0.6344340 

0.6335298 

-0.9041950E— 03 

4504400 

0.5619733 

0.5612927 

-0.6805807E— 03 

5132720 

0.4 77 2685 

0.4769768 

-0.2917275E-03 

576 1040 

0.3803502 

0.3804867 

0.1364686E-03 

6389360 

0.2713221 

0.2717314 

0.4093461E-03 

7017670 

0.1503622 

0.1507260 

0.3638361E-03 

7645990 

0. 1772289E-01 

0. 1786852E-01 

0.1456330E— 03 

8274310 

-0.1262693 

-0.1263368 

-0.6754883E-04 

8902630 

-0.2812140 

-0.2813877 

-0.1736917E-03 

9530950 

-0.4466380 

-0.4467816 

-0. 14362486-03 

0159260 

-0.6219965 

-0.6220134 

-0.1694 262 E— 04 

0787580 

-0. 8^66750 

-0.8065863 

0.88 75132 E-04 

141 5900 

-1.0000000 

-0 

1.0000000 
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TABLE n. - COMPARISON OF y*' WITH f'(x) 


X 

dy/dx 

dy/dx 

dy*/dx - dy/dx 

0 

0 

-0 . 6664473E-02 

-0.6664473E-02 

0. 6283180E-01 

0.1254983 

0. 1261036 

0.6053094E-03 

0. 1256636 

C. 2500058 

0.2538756 

0.3 869 805E-02 

0. 1886-954 

0.3725378 

0.3766517 

0.4 1 13 883E-02 

0.2513272 

0.4921210 

0.4944317 

0.2310704E-02 

0. 3141590 

0 • 6C77997 

0.6072157 

-0. 5840361 E-03 

0.3769908 

0.7186415 

0.7150037 

-0.3637798E-02 

0.4398226 

0.8237424 

0.8177956 

-0.5946755E-02 

0. 5026544 

0.9222328 

0.9155916 

— 0.6641 194E-02 

0. 5654862 

1 .0 132823 

1.0083916 

-0.4890755E-02 

0.6283180 

1.0961049 

1.0961955 

0 . 9052455E -04 

0.6911498 

1.1699639 

1.1752873 

0.5 323440E-02 

0. 7539816 

1.2341759 

1.2419509 

0.7775053E-02 

0. 8168134 

1.2881158 

1.2961864 

0 . 8070603E-02 

0 . 8796451 

1.3312202 

1.3379936 

0. 6 77 3 457 E -02 

0.9424770 

1.3629912 

1.3673728 

0.4381 567E— 02 

1.0053090 

1.3829996 

1.3843237 

0.1 324087E-02 

1.068 1410 

1.3908878 

1.3888464 

-0.2041385E-02 

1. 1309720 

1.3863723 

1.3813069 

-0.5065411E-02 

1.1938040 

1. 3692459 

1.3620709 

-0.7174999E— 02 

1.2566360 

1.3393793 

1.3311386 

-0.8240774E-02 

1.3194680 

1.2967224 

1.2885097 

-0 .8212730E-02 

1.3823000 

1.2413053 

1.2341844 

-0.7120892E-02 

1.4451310 

1.1732398 

1.1681638 

-0.5075917E-02 

1. 5079630 

1.0927146 

1.0904458 

-0.2268776E-02 

1.5707950 

1.0C00021 

1.0010314 

0.1029 342 E-02 

1.6336270 

0.8954524 

0.8999205 

0 . 446809 1 E-02 

1.6964590 

0.7 794939 

0.7871132 

0.7619314E— 02 

1. 7592900 

0.6526328 

0.6626115 

0 .99 78689 E—0 2 

1.8221220 

0.5154434 

0.5264115 

0.1 096 810E-0 1 

1. 8849540 

0 .3685770 

0.3785150 

0.9937912E— 02 

1.9477860 

3.2127518 

0.2189220 

0.61 70245 E-02 

2.C106180 

0.4875046E-01 

0.4763266E-01 

— 0.1117799E— 02 

2.0734490 

-0.1225800 

-0.1305877 

-0.8007 71 3E-02 

2. 1362810 

-0. 3C03440 

-0.3109823 

-0. 1063830E-01 

2. 199 1130 

-0.4835848 

-0.4935481 

— 0.9963315E— 02 

2.2619450 

-0.6713007 

-0.6782853 

-0.6984554 E-02 

2. 3247770 

-0.8624470 

-0.8651937 

-0.2746679E-0 2 

2.3876080 

-1.0559385 

-1.0542704 

0. 1668 110E-02 

2. 4504400 

-1.2506672 

-1.2455219 

0.5145237E-02 

2. 5132720 

-1.4454896 

-1.4385689 

0.6920710E-02 

2.5761040 

-1.6392453 

-1.6330362 

0.6209150E-02 

2.6389360 

-1.8307597 

-1. 8289229 

0.1836717E-02 

2. 701 7670 

-2.0188470 

-2. 0215197 

-0.2672 642E-02 

2.7645990 

-2.2023296 

-2. 2061239 

-0.3 794 28 3 E-02 

2.8274310 

-2. 3800253 

-2.3827324 

-0.2707154E-02 

2.8902630 

-2.5507663 

-2.5513453 

— 0.5789 101E-03 

2.9530950 

-2.7134030 

-2.7119624 

0.1 440 585 E -02 

3.0159260 

-2.8668072 

-2. 8645815 

0.2225727E-02 

3.0787580 

-3.0C98889 

-3.0092074 

0.6815 19 5E — 03 

3. 141 5900 

-3.1415873 

-0 

3.1415873 
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TABLE m. 

- COMPARISON OF 

/" x f 

/ 1 y* dx WITH 

x o 

f Xf f (x) dx 
x o 

X 

/y dx 

fy* dx 

fy dx - fy* dx 

0 

-0 

0.4712424E-09 

0.4712424E-09 

0.6283180E-0I 

-0.6274915E-01 

-0. 6 27 5994 E— 01 

-0.1079123E-04 

0. 1256636 

-0. 1250032 

-0.1250237 

-0.2050772E-04 

0. 1884954 

-0.1862709 

-0.1862868 

-0.1594424E-04 

0.2513272 

-0.2460688 

-0.2460646 

0.4187 226E-05 

0.3141590 

-0 .30392 52 

-0.3038922 

0.3308058E-04 

0.3769908 

-0 .3593838 

-0.3593241 

0.5961955E-04 

0.4398226 

-0.4120071 

-0.4119351 

0.7204339E-04 

0.3026544 

-0.4613806 

-0.4613191 

0.6152 689E-04 

0.5654862 

—0 .507 1 1 57 

-0.5070902 

0.2558529E-04 

0.6283180 

-0.5488533 

-0.5488819 

-0.2858788E-04 

0.6911498 

-0. 5862666 

-0.5863489 

-0.8233 63 7E-04 

0. 7539816 

-0.6190642 

-0.6191802 

-D.1159683E— 03 

0. 8168134 

-0.6469929 

-0.6471125 

-0.1196191E-03 

0. 8796451 

-0.6698399 

-0.6699318 

-0.9193271E-04 

0.9424770 

-0.6874352 

-0.6874730 

-0.3786385F-04 

1.0053090 

-0.6996535 

-0.6996202 

0.332891 9E-04 

1.0681410 

— 0.7C64159 

-0. 7063063 

0.1 095 608 E-0 3 

1.1309720 

-0.7C76914 

-0.7075135 

0.1778 826 E-03 

1.1938040 

-0.7034979 

-0.7032713 

0.2265275E-03 

1.2566360 

-0.6939030 

-0.6936558 

0.2471 805E-03 

1.3194680 

-0.6790246 

-0. 6787889 

0 .2356768E-03 

1. 3823000 

-0.6590311 

-0.6588390 

0.1920 834E-0 3 

1.4451310 

-0.6341418 

-0.6340211 

0. 1206994E-03 

1.5079630 

-0.6046244 

-0.6045948 

0.2 952 66 5 E-04 

1.5707950 

-0.5707971 

-0.5708675 

-0.7043779E-04 

1.6336270 

-0.5 330258 

-0.5331921 

-0.1662821E-03 

1.6964590 

-0.4917232 

-0.4919678 

-0.2445988E-03 

1.7592900 

-0.4473476 

-0.4476407 

-0.2930835E-03 

1.8221220 

-0.4 CO 3982 

-0.4007008 

-0.3026240E-03 

1. 8849540 

-0.3514170 

-0.3516866 

-0.2695 32 2E-03 

1.9477860 

-0.3009838 

-0.3011819 

-0.1981072E-03 

2.0106180 

-0.2497133 

-0.2496167 

-0.1034811E-03 

2.0734490 

-C. 1982536 

-0.1982667 

-0.1311488E-04 

2. 1362810 

-0.1472790 

-0. 1472320 

0.4703179E-04 

2. 1991130 

— 0.9749204E-01 

-0.9742580E-01 

0.662431 1 E-04 

2.2619450 

-9.4961565E-01 

-0 . 495687 6E-01 

0.4689349E-04 

2. 3247770 

-0.4390597E-02 

-0.4390209E-02 

0.3881287E-06 

2.3876080 

0.3742823E-01 

0.37371 32 E — 01 

-0.5691033E-04 

2. 4504400 

0.7507864E-01 

0.7497072E— 01 

-0.1D79114E-03 

2.5132720 

0 .1077915 

0. 1076525 

-0.1390344E— 03 

2.5761040 

0.1 347983 

0.1346546 

-0. 1437012E-03 

2.6389360 

0.1553342 

0.1552091 

-0.1250263E-03 

2. 701 7670 

0.1686435 

0. 1685442 

-0.9 92 5 66 3 E-0 4 

2. 7645990 

0. 1739844 

0. 1739015 

-0.8292 124E-04 

2.8274310 

0.1706328 

0.1705520 

-0.8086115E— 04 

2.8902630 

0. 1578875 

0. 1577984 

-0.8911267E-04 

2.9530950 

0.1350748 

0. 1349750 

-0.9971 48 5 E-04 

3.0159260 

0.1015535 

0. 1014485 

— 0.1050094E-03 

3.0787580 

0.56717 49E-0 1 

0.566152 1 E— 0 1 

-0.1022718E-03 

3. 1415900 

0.2652407F— 05 

-0.9539165F-04 

-0.9804405E— 04 
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